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13  ABSTRACT 


The  objective  of  this  research  program  was  development  of  finite  element  pro¬ 
cedures  for  prediction  of  stresses  and  deformations  in  the  vicinity  of  underground  excavation. 

Two  models  of  rock  behavior  we  re  selected.  In  one  the  rock  is  treated  as  isotropic 
elastic-plastic  following  a  generalized  Mohr-Coulomb  law  and  in  the  other  the  rock  is  iso¬ 
tropic  elastic-brittle  subject  to  Griffith  or  modified  failure  theory. 

For  each  model,  mathematical  relationships  governing  stress -strain  behavior  and  pro¬ 
gressive  failure  were  developed.  Finite  element  computer  programs  incorporating  each  of  the 
two  models  were  coded.  Preliminary  to  this  development,  a  revised  version  of  Zienkiewicz's 
no-tension  analysis  was  programmed. 

The  procedures  developed  all  for  initial  stresses  in  rock,  arbitrary  shape  and  size  of  the 
opening,  any  given  sequence  of  construction/excavation,  material  nonhomogeneity  and  pro¬ 
gressive  failure. 

This  report  is  in  three  parts:  Volume  1-Mhin  Document;  Volume  2-Computer  Program 
tier's  Manual;  Volume  3-Computer  Programs. 

Volume  1,  Main  Document,  contains  the  theoretical  development,  and  discussion  of  appli¬ 
cation  of  the  techniques  to  theoretical  solutions,  experimental  data.,  published  results  and  a 
parametric  study  of  forces  in  structural  supports  for  a  tunnel  excavation. 
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TECHNICAL  REPORT  SUMMARY 


Program  Objectives 

The  objective  of  this  research  program  was  development  of  finite  element 
procedures  to  prediet  stresses,  deformations  and  progressive  failure  of  roek  associ¬ 
ated  with  underground  excavations.  For  applicability  to  arbitrary  sequence  of  exca¬ 
vation  operations,  it  was  neeessary  that  the  procedures  developed  allow  for  arbitrary 
initial  stresses  in  roek,  arbitrary  size  and  shape  of  the  opening  and  progressive  fail¬ 
ure.  Plane  strain  conditions  and  two  different  types  of  material  behavior  were  con¬ 
sidered.  Roek  was  treated  as  an  isotropie  elastic-plastie  generalized  Mohr-Coulovnb 
material  in  one  model  and  as  an  elastic-brittle  material  following  Griffith  theory  of 
fraeture  in  the  other. 

Background 

In  previous  applications  of  the  finite  element  method  to  roek  mechanics,  elastie- 
plastie  behavior  of  roek  has  been  modeled  as  nonlinear  elastic  for  computational  con¬ 
venience.  Further,  it  was  assumed  that  the  results  of  a  one-dimensional  test  could 
be  generalized  to  three-dimensional  analysis  t  i rough  the  use  of  an  equivalent  stress- 
equivalent  strain  curve.  In  some  applications,  two  stress  or  strain  parameters  were 
used.  These  procedures  are  unsatisfactory.  Assumption  of  isotropie  elasticity  assumes 
that  the  principal  directions  of  stress  and  strain  coincide.  In  plastieity  this  is  not  true. 
Also,  roek  behavior  is  characterized  by  a  significant  part  of  deformation  being  irre¬ 
versible.  For  this  reason,  the  meehanieal  behavior  in  unloading  is  different  from  that 
in  loading.  For  roek  with  preexisting  joints  or  developing  tensile  eracks,  a  'no  tension' 
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procedure  is  often  adopted.  In  this  method,  a  ik-^ar  elastic  solution  is  obtained  and 
all  tensile  stress  redistributed  simultaneously.  Actually,  as  cracking  progresses, 
the  rock  on  either  side  of  the  crack  is  relieved  of  stress  and  a  stress  concentration 
develops  near  the  crack  tip.  Conventional  procedures  ignore  these  effects  and  the 
progressive  nature  of  crack  development,  leading  to  erroneous  conclusions  regarding 
stresses  around  underground  openings. 

Accomplishments  Under  the  Present  Program 

The  research  conducted  under  this  contract  has  resulted  in  development  of 
computer  programs  based  on  more  realistic  simulation  of  material  behavior.  Thf 
incremental  theory  of  plasticity  has  been  used  to  characterize  the  stress-strain  be¬ 
havior  of  elastic-plastic  rock.  Role  of  kinematic  constraint  of  plane  strain  in  develop¬ 
ment  of  residual  stresses  in  rock  has  been  examined  on  the  basis  of  Hill’s  theory. 

New  techniques  have  been  developed  for  study  of  initiation  and  propagation  of  fracture 
in  rock  following  Griffith’s  theory  or  the  modified  Griffith  theory.  Allowing  for  sequen¬ 
tial  fracture  of  various  elements  in  a  system,  the  effect  of  progressive  stress  redis¬ 
tribution  in  the  remaining  system  is  correctly  incorporated.  Arbitrary  initial  stress 
states,  arbitrary  sequence  of  excavation  (or  construction),  arbitrary  size  and  shape 
of  opening,  and  nonhomogeneous  material  properties  were  allowed  for.  The  actual 
construction  operations  can  be  simulated.  The  procedures  developed  were  applies  to 
several  typical  problems  in  rock  mechanics  as  well  as  to  some  theoretical  and  labora¬ 
tory  studies  for  the  purpose  of  verification  and  illustration.  These  were  used  to  carry 
out  parametric  studies  to  examine  the  influence  of  rock  oropertie  upon  th ;  stresses 
in  steel  supports  in  a  tunnel. 
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POD  Implications 


The  procedures  developed  provide  useful  means  for  study  of  stability  of  under¬ 
ground  excavations  based  on  stresses  and  deformations  associated  with  the  mining 
operations,  structural  support  evaluation,  safety  analyses  of  openings,  study  of 
blasting  effectiveness  under  certain  conditions,  evaluation  of  mining  sequences,  study 
of  vulnerability  and  serviceability  of  underground  structures  etc. 

Organization  of  the  Report 

This  report  is  in  three  parts  as  follows: 

Volume  1  -  Main  Document 

Volume  2  -  Computer  Program  User's  Manual 

Volume  3  -  Computer  Programs 

Volume  1  contains  the  main  body  of  the  report  including  the  theoretical  development,  pro¬ 
gram  verification  and  case  studies.  Chapter  I  reviews  previous  efforts  in  the  general 
research  area  and  describes  the  objectives  and  methods  of  the  present,  research  in  the 
historical  context.  Chapter  II  describes  the  mechanical  behavior  of  rock  and  the  ideali¬ 
zations  used  in  the  research  under  report.  The  basis  and  methods  of  the  finite  element 
theory  are  briefly  discussed  in  Chapter  III  leading  to  the  formulation  of  matix  equations,, 
Chapter  IV  gives  details  of  the  analysis  technique  for  isotropic  elastic-plastic  generalised 
Mohr-Coulomb  rock  materials  and  Chapter  V  gives  the  numerical  analysis  procedure  for 
jointed  rock  and  rock  subjected  to  progressive  fracture  following  Griffith  or  modified 
Griffith  theory.  Examples  of  application  are  included  in  Chapters  IV  and  V.  Chapter  VI 
presents  application  of  the  elastic-plastic  analysis  computer  program  to  a  parametric 
study  to  evaluate  the  influence  of  rock  properties  on  stresses  in  steel  supports  for  specified 
initial  stresses  and  design  of  the  opening. 
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In  the  original  proposal,  model  testing  to  verify  some  aspects  of  rock  behavior 
under  plane  strain  conditions  was  foreseen.  The  effort  under  the  present  contract 
covered  procurement  of  suitable  plane  strain  test  equipment  and  design  of  suitable 
test  material.  Appendix  B  includes  a  report  on  this  effort. 

Volume  2  of  the  report  contains  description  of  the  three  computer  programs 
developed  under  the  contract  along  with  fortran  listings  and  instructions  for  input  pre¬ 
paration.  The  input  definition  and  the  listings  are  for  the  IBM  370/165  version. 

The  programs  are  the  primary  content  of  volume  3.  These  are  available  on 
magnetic  tape  from  DDC-TC,  U.S.  Department  of  Commerce,  Springfield,  Virginia 
22151,  telephone  (703)  321-8517. 
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CHAPTER  I.  INTRODUCTION 


Tne  terrestrial  crust  is  in  a  complex  state  of  stress.  Excavations  in  this 
stressed  medium  profoundly  influence  the  distribution  of  stress  which  in  turn  de¬ 
termines  the  stability  of  the  pit,  slope  or  underground  opening  as  the  ease  may  be. 
Feasibility  ol  a  project  is  related  to  economics  of  safe  and  stable  construction. 

As  an  aid  to  decision  making,  it  is  necessary  to  develop  quantitative  information 
regarding  the  'state'  of  insitu  rook  and  predict  the  stability  characteristics  of  rock 
under  a  change  in  mechanical  environment  associated  with  the  proposed  construction 
or  excavation  operations.  Accordingly,  engineeiing  investigations  in  rock  mechan¬ 
ics  are  motivated  by  the  following  two  objectives: 

i.  Evaluation  of  the  'state'  of  insitu  rock.  This  includes, 
among  others,  investigation  of  pre-existing  stresses 
in  the  rock  and  mapping  of  discontinuities,  material 
symmetries  and  other  physical  data  relevant  to 
engineering  decision  process. 

ii.  Evaluation  of  changes  from  the  initial  state,  asso¬ 
ciated  with  excavation  or  construction  operations, 
including  among  other  factors,  changes  in  the  stress 
field,  along  with  consequent  yield  of  material,  ex¬ 
tension  of  existing  faults  and  discontinuities  and 
development  of  new  ones. 

The  research  effort  covered  by  this  report  was  directed  towards  the  second 

objective.  In  order  to  realistically  predict  the  stresses,  deformations  and  distress 

\ 

of  rock,  it  is  necessary  that  the  method  of  analysis  used  take  account  of  the  following- 
i.  The  initial  state  including  pre-existing  stresses, 
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discontinuities,  nonhomogeneity,  and  material 
symmetries,  if  any. 

ii.  The  sequence  of  excavation  or  construction  operations. 

lii.  The  mechanical  behavior  of  rock  i.e. ,  deforma¬ 
tion,  yield,  failure  or  fracture  of  rock  under 
changes  in  stress  environment. 

iv.  Interaction  with  structural  supports,  if  any,  used 
to  improve  safety  and  stability. 

Traditionally,  limit  equilibrium  methods  have  been  used  to  predict  stress  states 
and  factors  of  safety.  The  theory  of  elasticity  has  been  extensively  used  to  evalu¬ 
ate  stresses  in  rock  in  the  vicinity  of  underground  openings.  However,  these  meth¬ 
ods  are  applicable  only  for  linear  elasticity,  homogeneous  materials,  simple  geome¬ 
try  (e.g.  circular  or  elliptical  openings),  and  isotropic  materials.  None  of  these 
assumptions  are  valid  for  rock.  Furthermore,  the  excavation  was  assumed  to  be 
carried  out  in  a  single  step.  Influence  of  sequential  nature  of  any  construction 
could  not  be  taken  into  account  by  these  methods. 

With  the  advent  of  high  speed  digital  computers,  numerical  methods  of  analy¬ 
sis  have  proved  to  be  increasingly  useful.  Of  available  techniques,  the  most  power¬ 
ful  is  the  finite  element  method.  The  method  was  introduced  by  Clough  (1960)  and 
was  essentially  an  offshoot  of  the  general  matrix  structural  analysis  techniques 
developed  in  the  aerospace  indusiry  (Argyris  (1960),  Turner  et  al.  (1956)  ).  Develop¬ 
ment  f  the  method  has  been  rapid  and  a  considerable  volume  of  literature  has  accu¬ 
mulated  over  a  relatively  short  period.  Among  the  more  important  references  are 
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the  text  books  by  Zienkiewicz  (1972),  Glen  (1972),  Desai  (1972)  and  proceedings 
of  conferences  and  symposia  (WPAFB  (3965,  1968,  1971),  Vanderbilt  (1969),  U.S.- 
Japan  Seminar  (1969),  IUTAM  Symposium  (1970)  ).  Mathematical  foundations  of  the 
method  (Zlamal  (1968),  Oliviera  (1968),  Oden  (1969),  Aubin  (1972)),  its  relationship 
to  variational  methods  (Melosh  (1963),  Pian  and  Tong  (1968)),  convergence  of  se¬ 
quences  of  approximate  solutions  (Walz  et  al.  (1968),  Yamamoto  and  Tokuda  (1971), 
Key  (1966),  Oliviera  (1968))  have  received  attention. 

The  finite  element  method  has  been  applied  to  rock  mechanics  problems 
and  its  special  suitability  stems  from  the  fact  that  any  geometrical  configuration 
can  be  considered  and  nonhomogeneity  of  material  does  not  present  any  difficulty 
in  the  solution  process.  Also  arbitrary  loads,  including  body  forces  and  surface 
loads,  and  displacements  can  be  prescribed  for  the  problem. 

King  (1965)  developed  finite  element  techniques  for  incremental  construction 
which  were  used  by  Clough  and  Woodward  (1965)  and  Sandhu  et  al.  (1967).  King 
also  incorporated  a  technique  for  relieving  stresses  in  individual  elements  in  analy¬ 
sis  of  time -dependent  problems.  This  procedure  has  been  since  used  to  simulate 
excavations  by  Nair  et  al.  (1968),  Dunlap  and  Duncan  (1969),  Duncan  and  Chang  (1970), 
among  others.  It  also  forms  the  basic  proc  edure  in  analysis  of  rock  as  'no  tension' 
material  by  Zienkiewicz  et  al.  (1968)  and  also  in  Zienkiewicz  et  al.'s  (1968)  ana'ysis 
of  elastic-plastic  materials.  Kulhawy  (1972)  improved  the  procedure  somewhat 
by  using  stresses  at  the  excavated  surface  rather  than  stresses  at  centers  of  surface 
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elements  to  simulate  stress  relief  due  to  excavation. 

Nonlinear  material  behavior  was  treated  by  Gallagher  ct  al.  (1962)  using 
the  'initial  strain'  approach.  Finite  element  analysis  procedures  were  developed 
for  nonlinear  elasticity  by  several  investigators.  The  approach  v. as  to  approximate 
a  one-dimensional  stress-strain  curve  by  a  mathematical  relationship.  Wilson  (1965) 
used  the  bilinear  law.  Multi -linear  laws  have  been  proposed  by  Dunlap  et  al.  (1968). 
Salmon,  Berkc  and  Sandhu  (1969)  incorporated  Richard-Goldberg  (1965)  and  Ram- 
berg-Osgood  law  (1943)  and  Wilson's  bilinear  representation  in  a  single  program. 
Kulhawy  and  Duncan  (1969)  and  Kulhawy  (1972)  adopted  Kondner's  (1963)  hyperbolic 
representation  of  the  stress-strain  curve.  These  mathematical  equations  use  one 
or  more  parameters  to  get  best  fits  and  are  unsatisfactory  insomuch  as  the  slope  of 
the  curve  is  not  closely  approximated.  Desai  (1971)  used  spline  functions  to  obtain 
the  st  fit.  This  way  he  was  able  to  obtain  the  best  fit, not  only  to  the  stress-strain 
plot, but  also  to  the  slope  of  the  curve  which  is  the  information  actually  used  in  the 
computer  program.  Derai  (1972)  extended  his  work  to  allow  for  the  effect  of 
confining  pressure.  This  was  achieved  by  using  bi-eubie  splines.  Singh  and  Chang 
(1972)  further  developed  upon  Dcsai's  work  to  use  spline  fit  by  approximating  curves 
independent  of  stress  path  and  calculated  secant  modulus  and  poisson's  ratio.  The 
procedure  has  been  extended  to  laminates. 

To  extend  the  procedure  to  three-dimensional  analysis,  the  general  approach 
has  been  to  interpret  the  one-dimensional  test  data  as  equivalent  stress-equivalent 
strain  plot.  Second  invariant  of  the  stress  derivative  and  the  octahedral  shearing 
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strain  are  the  quantities  most  often  used.  Defining  a  shear  modulus  in  this  manner, 


the  Poisson's  ratio  is  assumed  constant.  Kulhawy  (1972)  defined  both  elastic  modu¬ 


lus  and  Poisson's  ratio  as  stress-dependent  quantities.  This  approach  assumes  iso¬ 


tropic  behavior  throughout.  Another  approach  using  the  laws  of  plasticity  was  used 


by  Salmon,  Berke  and  Sandhu  (1969)  and  Sutherland  (1970).  All  these  approaches 


are  essentially  empirical.  If  the  material  is  indeed  nonlinear  clastic,  one  can  either 


treat  it  as  Green-elastic  and  set  up  on  energy  functional  for  isothermal  or  adiabatic 


conditions  as  the  ~ase  might  be,  or  directly  treat  the  stress  tensor  as  a  nonlinear 


function  of  the  strain  tensor  and  for  isotropic  materials  using  the  Cayley-Hamilton 


theorem  to  obtain  a  representation.  Evans  and  Pister  (1966)  showed  that  for  energy 


functional  cubic  in  strain,  five  constants  are  required  to  define  the  stress-strain 


law,  and  these  cannot  be  obtained  from  any  single  test.  Similarly  the  constants  in 


any  Cayley-Hamilton  representation  cannot  be  obtained  from  a  single  test.  It  is 


well-known  that  equivalent  stress-equivalent  strain  plots  obtained  from  different 


tests  on  rocks  and  soils  are  path  dependent  and  therefore  different. 


For  finite  element  analysis,  elastic-plastic  materials  (rate  independent 


materials  exhibiting  path  dependent  behavior  above  certain  stress  level)  were  treated 


as  nonlinear  elastic  by  several  workers.  In  rock  mechanics  the  deformation  theory 


of  plasticity  has  been  used  (Jaeger  and  Cook  (1969),  Brady  (1970),  Malina  (1970)). 


This  theory  assumes  the  principal  strain  directions  to  coincide  with  the  principal 


stress  directions,  assumes  no  volume  change  during  plasticity,  and  is  inapplicable 
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to  nonmonotonic  loading.  Hill  (1950)  has  discussed  the  shortcomings  of  this 

theory  in  detail.  It  suffices  to  note  here  that  this  theory  is  unsound  and  for 
this  reason  it  was  not  considered  in  the  present  research  program  as  a  possible 
model  for  mechanical  behavior  of  rock. 

Incremental  theory  of  plasticity  was  used  by  Zienkiewicz  et  al.  (1968), 
Swedlowet  al.(l965),  Marcal  (1969),  Marcal  and  King  (1967).  Drueker's  method 
for  evaluation  of  A  in  the  -normality  rule'  was  used  leading  to  the  tangent  modulus 
approach  for  von  Mises  materials.  Felippa  (1966)  developed  stress-strain  relations 
in  terms  of  total  incremental  strain  (elastic  and  plastic)  and  incremental  stress. 
This  relation  is  identical  to  the  one  proposed  by  Naghdi  (1960)  and  Prager  (1949) 
for  work-hardening  materials  but  lias  the  additional  advantage  that  it  is  valid  for 
perfectly  plastic  solids  for  which  Naghdi-Prager  relationships  are  undefined. 

Reyes  (1965)  and  Reyes  and  Deere  (1966)  used  a  rate  of  work  equation  to  develop 
stress-strain  relations  for  Mohr-Coulomb  materials.  This  approach  has  the  limi¬ 
tation  that  it  cannot  be  generalized  to  other  failure  laws.  Yamada  (1968)  used  an 
energy  rate  approach  for  von  Mises  materials.  These  energy  rate  approaches  are 
limited  in  application  and  can  be  derived  trom  Felippa's  more  general  formulation. 
Zienkiewicz  et.al.  (1968)  used  tdis  approach  in  conjunction  with  an  'initial  stress- 
technique  to  solve  for  stresses  without  interest  in  displacements.  Baker,  Sandhu 
and  Shieh  (1969)  used  Reyes  and  Deere's  formulation  for  finite  element  analysis 
oi  stresses  and  deformations  in  rock.  It  was  noticed  that  under  plane  strain 
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conditions,  a  radial  path  eannot  be  followed  in  the  Westergaard-Haigh  stress  spaee. 

Discontinuities  in  rock  appear  as  joints,  fissures  and  faults.  Modeling  of 
discontinuities  has  received  considerable  attention.  However,  exclusively,  the 
geometry  of  the  eraeks  and  joints  was  assumed  to  be  known.  With  this  as  the 
starting  point,  research  effort  has  principally  been  concentrated  on  eharaeteriza- 
tion  of  stress-strain  behavior  of  jointed  roek  and  on  incorporating  limiting  shear 
strength  of  joints  in  finite  element  analyses.  Anderson  and  Dodd  (1966)  considered 
open  joints  or  faults  as  surfaces  with  no  resistance  to  shear,  ineapable  of  with¬ 
standing  tensile  stresses  normal  to  their  planes.  This  capability  is  now  routinely 
incorporated  in  most  finite  element  programs.  A  two-dimensional  'soft'  material 
element  has  long  been  used  to  represent  weak  joint  planes  in  roek.  Duncan  and 
Goodman  (1968)  objeet  to  this  on  the  basis  of  the  large  number  of  elements  needed 
to  ensure  a  reasonable  'aspect  ratio'  in  the  shape  of  elements.  This  becomes 
a  problem  for  elements  representing  very  thin  joints.  For  rock  systems  having 
thin,  elosely  spaeed  parallel  joint  planes,  an  equivalent  orthotropic  elastie  eon- 
tinuum  representation  regarding  the  orientation  of  the  joints  as  the  plane  for  reflec¬ 
tive  symmetry  in  the  material  was  proposed  by  Dunean  and  Goodman  (1968).  A 
one-dimensional  element  with  shear  and  normal  stiffness  eharaeteristies  was 
developed  by  Goodman  et  al.  (1968).  Recently,  Heuze  et  al.  (1971)  have  intro¬ 
duced  nonlinear  meehanieal  properties  in  this  element.  Christian  is  credited  with 
development  of  an  element  capable  of  simulating  constant  shear  and  residual  shear 


characteristics.  All  these  investigations  were  eoneemed  with  setting  up  stress- 
strain  relationships  for  rock  with  pre-existing  joints.  Propagation  of  fracture 

has  been  considered  by  several  investigators.  Chan  et  al.  (1970) ,  Gross  et  al. 

(1968)  sought  to  calculate  stress  intensity  factors  at  erack  tips  to  determine 
whodier  a  given  eraek  would  propagate.  To  improve  aceuraey,  Wilson  (1971), 

Byskov  (1970)  and  Levy  (1971)  introduced  stress  singularity  elements  at  eraek  tips. 
Plan  (1971)  used  hybrid  finite  elements  to  obtain  improved  stress  intensity  factors 
without  recourse  to  singularity  elements.  Throughout,  it  was  assumed  the  geo¬ 
metry  of  the  eraek  was  known  beforehand.  These  procedures  eannot  be  used  to 
study  discontinuities  arising  as  a  result  of  fracture  under  changes  in  stress  environ¬ 
ment.  Also  these  eannot  consider  propagation  of  fraeture.  To  use  the  same  pro¬ 
cedure  for  pre-existing  as  well  as  post-failure  eraeks,  it  is  necessary  to  allow 
cracks  and  joints  within  elements.  Then,  the  mesh  layout  is  more  flexible  and 
arbitrary  failure  laws  can  be  used.  Malina  (1970)  used  this  approach  to  study  failure 
along  joint  planes  and  then  went  on  to  compute  the  amount  of  slip  and  the  accompanying 
stress  redistribution  on  the  basis  of  deformation  or  slip  theory  of  plasticity.  Re¬ 
cently,  Bock  (1971)  has  used  Malina's  technique  to  study  propagation  of  faulting  in 
roek. 

Simulation  of  rock  behavior  after  cracking  develops  can  be  achieved  using  a  bi- 
modular  elastic  representation  (Sandhu,  1969).  However,  recognizing  the  eraek  to  be  a 
plane  of  material  symmetry,  the  material  is  bimodular orthotropic  withno  resistance 
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to  tension  in  the  direction  normal  to  crack.  For  cracks  in  two  independent  directions, 
an  element  under  uniform  stress  can  be  regarded  as  bimodular  isotropic.  The  'no 
tension'  method  used  by  Zienkicwicz  ct  al.  (1968)  is  an  iterative  procedure  of  obtain¬ 
ing  tension  free  stress  states.  The  philosophy  rcfleets  an  iterative  solution  of  a 
bimodular  orthotropie  problem.  However,  the  mechanics  of  solution  eliminates  only 
the  tensile  principal  stresses.  This  amounts  to  using  a  non-symmetrie  stress-strain 
relationship  and  is  open  to  objection.  This  situation  was  corrected  by  Sandhu  (1971). 

New  procedures  for  finite  element  analysis  of  elastic-plastie  Mohr-Coulomb 
materials  and  elastic-brittle  materials  following  Griffith's  theory  under  plane  strain 
conditions  were  developed  in  the  eourse  of  the  present  research.  Arbitrary  geo¬ 
metric  configurations,  complex  boundary  conditions,  nonhomogencity  of  roek  are 
allowed  for  as  is  usual  in  finite  element  techniques.  Variational  formulations  of 
the  field  equations  shows  that  arbitrary  initial  sti’esscs  and  pore-water  pressures, 
if  any,  can  be  directly  included  in  the  problem.  Stress-strain  relations  in  plas¬ 
ticity  follow  Reyes  and  Deere  (19GG)  and  Felippa's  (1966)  formulation.  At  early 
stages  of  work  (Sandhu,  1971),  it  was  assumed  that  in  plane  strain  plasticity,  the  clastic 
and  plastic  components  of  strain  on  planes  normal  to  the  axis  separately  vanish.  Thus 
a  sudden  jump  in  axial  stress  was  introduced  to  ensure  continuity  of  deformation. 
Later,  Hill's  (1950)  analysis  of  plane  strain  was  adopted  using  the  vanishing  of  total 
strains.  In  an  arbitrary  construction  sequence,  a  given  element  may  be  subjected 
to  excursions  to,  from  and  within  the  yield  surfaee  in  the  stress  spaee.  Capability 
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to  allow  for  loading  followed  by  unloading  has  been  incorporated  in  the  analysis. 
For  jointed  rock  subject  to  progressive  fracture,  Griffith  and  Modified  Griffith 
theories  were  used  to  predict  fracture  initiation.  Bimodular  clastic  orthotropy 
describes  the  post-fracturc  behavior  of  rock  and  an  iterative-incremental  process 
is  used  to  define  the  sequential  fracture  of  elements  in  a  system  giving  the  extent 
of  various  fractures.  The  new  pr  ccdure  differs  significantly  from  the  'no  ten¬ 
sion'  concepts.  The  latter  define  a  region  where  the  clastic  solution  indicates 
tensile  stresses  and  the  solution  process  neutralizes  all  tensions  simultaneously. 
This  is  incorrect  and  does  not  allow  for  the  stress  redistribution,  in  its  neighbor¬ 
hood,  caused  by  the  formation  of  a  crack.  Sequential  cracking  constitutes  a  non¬ 
linearity  in  material  bshav1  or  and  linear  superposition  implied  in  simultaneous 
release  of  all  tensions  in  the  'no  tension'  procedure  is  not  valid. 

The  procedures  developed  allow  for  arbitrary  construction  or  excavation 
sequence,  and  can  allow  for  pre-existing  joints,  open  or  closed.  Linear  elastic 
elements  can  be  included  in  the  analysis.  The  procedures  can  include  tunnel  sup¬ 
ports  as  part  of  the  sequential  construction  scheme.  These  were  used  to  obtain 
numerical  solutions  to  several  problems.  Results  arc  included  in  the  report. 

Chapter  II  of  the  report  discusses  mathematical  modeling  of  mechanical 
behavior  of  rock.  In  Chapter  in  basis  of  the  finite  element  method  is  outlined. 
Chapters  IV  and  V  present  details  of  application  of  the  finite  element  method  res¬ 
pectively  to  plane  strain  analysis  of  elastic-plastic  Mohr-Coulomu  rock  and 
plane  strain  analysis  of  progressive  cracking  of  rock  following  Griffith's  theory. 
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Chapter  VI  reports  the  application  of  the  elastic-plastic  incremental  structure 
analysis  to  a  parametric  study.  The  objective  of  this  study  was  to  determine 
the  influence  of  rock  property  data  on  stresses  in  steel  supports  of  a  tunnel. 
Chapter  VII  contains  a  discussion  of  results  of  this  research  effort.  Appendices 
A  and  B  present,  respectively,  Griffith’s  theory  of  fracture  propagation  and  the 
experimental  research  efiort  under  the  contract. 
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CHAPTER  II.  MATHEMATICAL  MODELING  OF 
MECHANICAL  BEHAVIOR  OF  ROCK 


2.1.  Mechanieal  Behavior  of  Rock 

Figs.  II- 1  and  II— 2  show,  respectively,  typical  stress-strain  plots  for  a 
granite  and  a  marble.  Upon  loading  the  stress-strain  curve  is  almost  linear  and 
reversible  over  a  short  portion.  Unloading  from  higher  loads  docs  not  coincide 
with  initial  loading.  This  characteristic  along  with  rate  independence  distinguishes 
elastic-plastie  behavior.  Reloading  closely  follows  unloading  until  the  previous 
maximum  is  reached;  whereupon,  the  original  curve  is  followed.  This  leads  to 
some  simplifying  assumptions. 

i.  A  yield  point  exists  below  whieh  the  material 
is  linear  elastic. 

ii.  The  yield  point  corresponds  to  the  maximum 
stress  level  previously  attained. 

iii.  Unloading  and  reloading  paths  are  linear, 
coincident  and  parallel  to  the  initial  elastic 
loading  eurve. 

Fig.  II-3  shows  this  simplification.  Clearly  the  yield  point  ean  be  described  by 
the  permanent  or  irrecoverable  strain  or  the  area  bounded  by  the  loading  curve, 
the  unloading  curve  and  the  horizontal  axis.  Whereas  in  generalization  to  the  three- 
dimensional  ease,  the  stresses  and  strains  beeome  seeond  rank  tensors  and  are 
therefore  unordered;  the  area  is  still  a  scalar  product  and  retains  its  ordering 
characteristics.  To  this  extent,  it  is  often  preferred  as  a  measure  of  the  elastic 
limit.  Other  approaches  treat  the  elastic  limit  as  a  scalar  function  of  the  strain 


tensor. 


Stress-Strain  Curves  for  Cellar  City  Granite 
(Swanson,  1970) 


FIG.  H-2.  Stress-Strain  Curves  for 
(Swanson,  1!)70) 


Mechanical  behavior  of  rock  under  polyaxial  state  of  stress  has  been  examined 
in  the  light  of  brittle  failure  theories.  Four  regions  of  behavior  are  identified  in 
Fig.  II— 4 .  The  first  region  corresponds  to  closure  of  preexisting  open  cracks  and 
is  peculiar  to  compressional  loading.  In  region  II  material  behavior  is  linear  elastic. 
Fracture  initiation  occurs  near  the  end  of  this  region  in  accordance  with  Griffith  or 
Modified  Griffith  theory.  This  stage  also  corresponds  to  onset  of  nonlinearity  in  the 
relationship  of  stress  to  volumetric  strain  (Brace,  19G6).  Stable  fracture  propagation 
characterizes  region  III.  In  region  IV,  unstable  fracture  propagation  results  in  strength 
failure  and  rupture.  Differerces  in  loading  and  unloading  behavior  are  observed  (Walsh, 
1965). 

We  have,  thus,  two  general  approaches  to  the  characterization  of  stress-strain 
behavior  of  rock.  One  follows  the  theory  of  elastic-plastic  solids  without  consideration 
of  micro-mechanics  of  the  system.  The  other  uses  Griffith  theory  or  Modified  Griffith 
theory  to  rebate  deformation  and  failure  to  initiation  and  propagation  of  fracture.  It 
has  been  observed  (Swanson,  1970)  that  Mohr-Coulomb  failure  law  applies  for  mode¬ 
rate  values  of  confining  pressure  and  that  at  low  confining  pressures,  failure  is  by 
rupture.  Contrary  to  plastic  behavior,  strength  of  material  drops  to  almost  zero  in 
the  direction  normal  to  the  crack  if  rupture  theory  is  followed.  Figs.  II-5  and  II-6 
depict  typical  relationships  of  failure  strength  and  post-failure  behavior  in  relation 
to  confining  pressures.  It  is  reasonable  to  assume  that  the  material  is  linear-elastic 
upto  yield  or  rupture,  as  the  case  may  be,  and  that  the  post-failure  behavior  only  is 
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governed  by  the  theory  used  to  define  failure.  In  the  present  research  program,  both 
the  elastic-plastie  Mohr-Coulomb  failure  theory  and  the  Griffith  theory  have  been  used 
for  analysis  of  stresses,  deformation,  and  failure  of  rock. 

2.2.  Characterization  of  Elastic-  Plastic  Behavior 

Several  approaches  have  been  used  for  formulation  of  elastic-plastic  behavior 
of  rock.  One  approach  is  to  treat  rock  as  nonlinear  elastic  material.  It  is  argued 
that  this  assumption  is  a  reasonable  one  in  case  unloading  docs  not  occur  and  loading 
is  'proportional'  (or  radial  in  the  stress  space).  Under  such  circumstances,  the 
state  of  stress  is  described  completely  by  a  single  parameter.  If  stress-strain  data 
for  the  actual  stress  path  are  available,  these  can  be  directly  used  to  predict  defor¬ 
mation  response  to  applied  loads.  This  approach  is  attractive  because  of  its  apparent 
simplicity.  However,  stress  path  dependence  of  material  behavior  results  in  diffi¬ 
culties  in  correct  representation  of  behavior  under  arbitrary  stress  changes.  Attempts 
have  been  made  to  represent  strain  invariants  as  functions  of  stress  invariants 
and  vice  versa,  thereby  permitting  generalizations  to  arbitrary  stress  paths.  It  is 
customary  to  use  a  stress  or  strain  dependent  Young's  modulus  and  Poisson's  ratio. 
Wilson's  bilinear  law  (1965),  Goldberg-Richard  equation  (1965),  Ramberg-Osgood 
equation  (1943),  Prager's  formulation  (1938),  the  work  of  Kondner  (1963),  Duncan  and 
Chang  (1970),  Kulhawy  (1971),  Vallabhan  and  Reese  (1968),  Desai  (1971)  and  recent 
work  by  Singh  (1972),  are  all  based  on  this  philosophy. 

In  general,  there  are  two  distinct  approaches  to  characterization  of  nonlinear 
elastic  materials.  One  is  to  assume  the  existence  of  energy  functional  such  that  its 
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derivative  with  respect  to  any  strain  component  is  the  corresponding  component 
of  the  stress  tensor.  The  other  directly  writes  stress  as  a  tensor  function  of  strain. 


Assuming  the  exist snce  of  strain  energy  functional,  Rivlin  (1960)  proposed 


a  ip  iy 
1  l2  *3 


(II-l) 


a,0,y 


where  Ij  =  I2  =  2  f  ii  +  *ii  fjj  "  *ij  *ij  nnd  I3  =  determinant  (5ij  +  2  fjj)  -  1. 
Toupin  and  Bernstein  (1961)  proposed 


U  2  Eijkl  fij  fkl  +  6  Eijklmn  fij  fkl  (mn  +  **• 


(n-2) 


Neuber  (1969)  used 


U  =  u  (<»1#  <j>2 ,  <t>3) 


(n-3) 


where 


<b ,  =  A. . 
rl  i)  13 

^2  =  Aijkl  fij  fkl 

^3  =  Aijklmn  'ij  'kl  (mn 


(H-4) 


Neuber's  formulation  includes  Toupin  and  Bernstein's  as  a  specialization.  For  isotropy 


A|j  ,  Ajj^j ,  A|jjc|mn  are  taken  as  identity  tensors.  Evans  and  Pister  (1966'*  proposed, 


for  isotropic  elasticity 


U  =  U  (  <pv  cp2,  <p3) 


(n-5) 


as  a  cubic  or  biquadratic  expression  in  ( For  biquadratic  expansion 


A11  eh  2  ,  A21 
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requiring  nine  constants  to  describe  the  energy  functional.  Here  <^  =  f  ^ 

2  fij  V  ^3=  |  fij  fjk  fki  ■  It  is  important  to  note  that  ( . . - *»  0  does  not  imply 

linearity  in  the  limit.  Using  a  cubic  expansion,  there  arc  five  constants. 

Alternatively,  writing  stress  directly  as  a  polynomial  function  of 
strain,  Cayley-Hamilton  theorem  leads  to  the  following  relationship  for  isotropic 
materials 

Tij  =  a  5i)+  b  f(j+c  <(k  <kj  (H-7) 

where  a,b,c  are  functions  of  strain  invariants. 

It  is  clear  that  the  characterization  of  material  behavior  based  on  a  single  type 
of  test  cannot  be  expected  to  completely  represent  nonlinear  elasticity.  Even  if  it 
were  possible,  the  analysis  could  not  simulate  elastic-plastic  behavior  in  which  the 
plastic  strain  increments  are  normal  to  the  stress  increments  in  loading.  It  is  nec¬ 
essary  therefore  that  elastic-plastic  behavior  of  rock  be  modeled  using  the  theory 
of  elastic-plastic  continua. 

A  general  theory  of  elastic-plastic  continua  was  presented  by  Green  and  Naghdi 
(1965).  Stress-strain  relationships  in  plasticity  and  thermo-plasticity  have  been  re¬ 
viewed,  among  others  by  Drucker  (1950),  Naghdi  (1960),  Koitcr  (1953)  and  Hill  (1950). 
In  general,  theoretical  concepts  have  been  based  on  generalization  of  material  behavior 
in  a  one-dimensional  test. 


The  deformation  theory  of  plasticity  relates  increments  in  components  of  the 
strain  tensor  to  corresponding  components  of  the  stress  deviation  tensor.  This  theory 


is  unsound  theoretically  and  fails  to  properly  explain  physical  behavior  of  materials. 
It  has  been  used  to  represent  material  behavior  under  proportional  or  almost  pro¬ 
portional  loading  (Budiansky  (1959),  Havner  (1969),  Ilyushin  (1945)).  Hill  (1950)  has 
examined  the  shortcomings  of  the  theory  in  his  text.  These  will  not  be  discussed 
here,  and  we  shall  confine  our  attention  to  the  incremental  or  rate  type  theory  of 
plasticity. 

Using  Prandtl's  idealization,  shown  in  Figure  II-3,  it  is  assumed  that  in  one¬ 
dimensional  loading: 

1 .  There  is  a  yield  point  such  that  any  loading  beyond 

this  level  is  accompanied  by  permanent  deformation  and 
any  stress  changes  below  Oyfeld  have  linear  reversible 
response  independent  of  prior  deformation  history. 

2.  The  yield  point  is  a  function  of  history  of  prior  permanent 
deformation. 

3.  The  stress-strain  behavior  is  independent  of  sign  of  the 
stress. 

In  generalization  to  arbitrary  stress  states,  it  is  customary  to  introduce  a  yield 
surface  in  the  six-dimensional  (three-dimensional  for  isotropic  materials  and  nine¬ 
dimensional  for  non-symmetric  stress  tensor)  stress  space.  Allowing  for  history 
of  deformation  and  material  properties,  the  yield  condition  is: 


f  (o\.  ,  J\.  k.  ) 

ij  *  *  ij  ’  i  ' 


(II— 8) 


where  e".j  are  components  of  the  plastic  strain  tensor  (irreversible)  and  k|  are  mate¬ 
rial  parameters  depending  upon  history.  The  mapping  f  in  Eq.  II-8  defines  a  func¬ 
tional  on  the  linear  vector  space  V  spanned  by  <r..,  «  ,  k.  ,  and  orders  the  space 
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V  such  that  'loading'  and  'unloading'  can  be  defined  analogous  to  the  one-dimensional 
test  where  the  set  of  stress  states  is  an  interval  on  the  real  line  and  is  naturally 
ordered.  Thus  f  <  0  will  imply  elastic  states  and  plastic  deformation  is  possible 
only  for  f  -  0.  We  shall  denote  the  interior  of  the  surface  defined  by  f  =  0  as  D. 

Thus  f(u)  <  0  for  uGD.  Assuming  only  one  material  constant  k,  the  time  rate  of  f  is 
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where  X  is  a  positive  scalar  which,for  rate  independence, must  be  homogeneous  of 


order  one  in  <r...  Fcr  von  Mises  materials,  Equation  (11-13)  reduces  to  the  well- 


known  Prandtl-Reuss  relationships.  To  evaluate  X,  Drucker  defined  equivalent  stress 


<re  and  equivalent  plastic  strain  rate  *"e  as  second  invariants  of  the  stress  deviation 


tensor  and  of  the  plastic  strain  rate  tensor  respectively.  These  invariants  are  proportional 


to  the  I<2  norm  on  the  linear  vector  spaces  spanned  respectively  by  the  components  Sjj 
of  the  stress  deviation  tensor  and  the  components  ( jj  of  the  plastic  strain  rate  tensor. 


Thus 
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where  the  components  of  the  stress  deviation  tensor  are  related  to  the  stress  tensor 


by  the  relationship 
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Then 
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Writing  *n  ,  the  slope  of  the  <re  ,  *"e  curve  as  H,  we  obtain  upon  substitution  of 

(  e 


Equation  (11-18)  in  Equation  (11-13), 
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This  formulation  was  used  in  the  so  called  tangent-modulus  methods  e.g.  Swedlow  and 
Yang  (1965).  Clearly  this  analysis  satisfies  the  condition  f  =  0  and  the  normality  rule 
but  fails  to  satisfy  f  -  0  in  plastic  deformation. 

Hill  (1950)  used  the  normality  rule  assuming  X.  to  be  a  fourth  rank  tensor  linear 

•  # 
in  0-^1  and  introduced  a  plastic  potential.  Using  normality  as  well  as  the  condition  f  =  0, 

Prager  (1949)  obtained  for  f  =  f(cr.j,  f "••) 

f 
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Naghdi  (1960)  also  obtained  the  above  equation.  This  formulation  breaks  down  for 

elastic  perfectly  plastic  solids  where  f  is  independent  of  e"  .  Felippa  (1966)  obtained 

X.  in  terms  of  ijj,  the  total  strain  increment  as  follows. 

The  incremental  stress  is  related  to  the  incremental  clastic  strain  by  the  equation 

17 ij  =  Eijkl  <'kl  (n"22> 

where  are  components  of  the  elastic  strain  rate  tensor  and  E-j^  is  the  fourth  rank 
isothermal  elasticity  tensor  satisfying  the  symmetry  properties 

Eijkl  =  Ejikl  =  Eij  lk  =  Eklij  (II_23) 
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Writing  =  *ki  -  f",.  and  using  normality  rule, 


(11-24) 


°"ij  Eijkl  'kl 


E....  X.  -11 

lJkl  A  <r,. 


(11-25) 


Also  f  (cr„  ,  f"„)  -  0  implies 
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Using  normality  rule,  Equation  (11-26)  gives 
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Substituting  Equation  (11-25)  in  Equation  (11-17), 
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Hence 
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Substituting  Equation  (11-29)  in  Equation  (11-25) 
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This  formulation  relates  the  total  stress  increment  to  the  total  strain  increment  and  is 
valid  for  all  cases  including  perfect  plasticity.  For  work-hardening  plasticity,  it  can 
be  shown  to  be  identical  to  the  Prager-Naghdi  formulation  of  Equation  (11-21). 

Using  rate  of  work  equations,  it  is  possible  to  evaluate  X  in  terms  of  stress 
rates  for  von  Mises  or  generalized  Mohr-Coulomb  materials.  Yamada  (1968)  used 
shear  deformation  energy  rate  to  set  up  incremental  stress-strain  relationships  for 
von  Mises  materials.  Assuming  to  occur  in  f  as  the  second  invariant  only,  a 
tangent  modulus  approach  was  introduced.  Reyes  (1965)  developed  the  incremental 
stress-strain  relations  for  generalized  Mohr-Coulomb  materials  using  the  total  energy 
rate.  The  general  forms  of  these  equations  can  be  shown  to  be  specializations  of  the 
more  general  formulation  developed  by  Felippa. 

The  stress-strain  relations  expressed  by  Equations  (Ii-22)  and  (11-31)  for  the 
elastic  and  el  as  tic -pi  a  Stic  deformation  respectively  have  to  be  specialized  for  any 
kinematic  constraints  that  may  exist  in  the  case  of  plane  strain,  for  example 

^3=0;  i  =  1,2,3  (11-33) 

Thus  in  the  elastic  domain,  if  Ej^  is  invertible  such  that 

Cijkl  =  <Eijkl>  1  (II-34) 

the  kinematic  constraint  is  expressed  by  a  linear  relationship  between  components 


of  incremental  stress  tensor, 


*13  Ci3kl  °kl 


For  the  clastic-plastic  deformation,  writing 


Cijkl  [Eijkl  <pkm  pln  "^lmn^] 
the  kinematic  constraint  is 

*i3  =  Ci3kl  °kl  ”  0 


(11-35) 


(11-36) 


(11-37) 


It  is  important  to  note  that  C^i  depends  upon  cr jj  and  .  Consequently  the 
relationship  expressed  by  Equation  (11-37)  is  nonlinear  and  involves  <r.j,  as  well 
as  (r^j .  For  a  perfectly  plastic  solid,  Cjg^j  is  a  function  of  stress  and  elastic  pro¬ 
perties.  For  a  von  Mises  solid  the  problem  was  studied  by  Hill  (1950).  Figure  II— 7 
shows  a  plane  strain  loading  and  unloading  cycle  for  (r^  =  =  0.  Path  OA 

represents  elastic  loading  to  yield  at  A.  ConU  lued  loading  requires  the  stress  path  to 
follow  AD, the  trace  of  the  loading  surface  on  the  <r^j,  <Tgg  plane.  Loading  from  A  to  D 
is  accompanied  by  development  of  plastic  strains  e"n,  e"2z  and  a  residual  stress 
0-33.  Upon  elastic  unloading  from  D  to  E,  the  path  is  parallel  to  OA.  The  intercept 
OE  represents  the  residual  stress  <r 33.  This  stress  remains  in  the  material  upon 
removal  of  the  two-dimensional  load.  Accompanying  this  residual  stress  are  elastic 
as  well  as  plastic  strains  satisfying  the  kinematic  constraints.  The  above  theory  was 
used  to  develop  mathematical  model  representative  of  plane  strain  elastic-plastic  be¬ 
havior  of  Mohr-Coulomb  rock.  The  details  are  contained  in  Chapter  IV  of  this  report. 
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Residual  Stresses  Associated  with  Elastic- Plastic  Deformation  Under  Plane  Strain 


2.3.  Mechanical  Behavior  of  Brittle  Rock 


Stress-strain  behavior  of  rock  under  polyaxial  compression  has  received  con¬ 
siderable  attention  (Brace  (1904),  Cook  (1905),  Ricniawski  (1907),  Brady  (1970)  ). 

In  a  one-dimensional  test  the  stress-strain  curve  can  be  divided  into  four  regions. 
Figure  II— 4  taken  from  Brady  (1970)  shows  the  different  regions.  These  are  char¬ 
acterized  by  the  following: 

i.  Closure  of  initial  flaws  under  compression 

ii.  Microcrack  initiation 

iii.  Stable  fracture  propagation 

iv.  Unstable  crack  propaga,son  leading  to  rupture 

In  compression,  closure  of  preexisting  microcracks  is  essentially  completed 
over  a  small  increment  of  stress.  A  nonlinear  stress-strain  relation  has  been  ob¬ 
served  for  this  zone,  showing  the  high  deformability  due  to  open  cracks  decreasing 
as  the  cracks  close.  Friction  is  developed  when  surfaces  of  the  initial  flaws  contact 
each  other,  resulting  in  an  increase  in  the  observed  modulus  of  elasticity. 

When  the  applied  stress  is  increased  further,  the  mechanical  behavior  as  seen 
from  the  stress-strain  curve  is  essentially  linear  elastic,  and  the  modulus  of  elasticity 
is  constant.  However,  some  sliding  across  the  faces  of  closed  cracks  may  occur. 
Walsh  (1965)  demonstrated  that  the  loading-unloading  process  at  this  stage  shows 
hysteresis  in  the  stress-strain  curve  (Figure  II-8).  Some  of  the  frictional  strength 
contributed  by  the  roughness  on  the  crack  surfaces  is  overcome  by  sliding,  and  upon 
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REGION  I 


unloading,  the  crack  does  not  move  back  to  the  original  position  immediately.  Thus, 
even  in  the  linear  elastic  region  II,  rocks  may  exhibit  irreversibility  as  well  as 
hysteresis  (Figure  II-l). 

Sliding  along  each  initial  crack  produces  new  cracks  at  the  ends  of  the  initial 
crack.  Tests  on  quartzite  led  Bieniawski  (1967)  to  conclude  that  the  relative  dis¬ 
placement  of  crack  surfaces  is  the  primary  factor  which  influences  fracture  initiation. 
In  region  III,  the  microcracks  propagate  with  increasing  stress  level.  Stable  fracture 
propagation  is  a  function  of  stress  only  and  the  process  is  quasi-static.  Initiation  and 
propagation  of  cracks  is  reflected  in  the  aeparture  from  linearity  of  the  stress-strain 
curve  (Brace,  (1966)  ).  Crack  growth  starts  right  after  the  initiation  of  new  surfaces. 
Tests  in  glass  and  rock  (Hoek  and  Bieniawski,  (1965)  )  have  shown  that  stable  fracture 
propagation  follows  a  curved  path  leading  to  a  direction  parallel  to  the  major  compres¬ 
sive  stress,  and  the  length  of  propagation  is  related  to  the  ratio  of  major  principal 
stress  to  minor  principal  stress  and  related  to  the  initial  crack  as  well  (Paul,  et.al. 
(1967)).  Pattern  of  crack  array  can  also  influence  'die  fracture  propagation. 

Before  the  applied  stress  reaches  its  peak,  some  abrupt  changes  on  propagation 
behavior  would  happen.  Although  the  stress  is  still  increasing,  the  slope  of  the  stress- 
strain  curve  is  decreasing.  As  soon  as  slope  reaches  its  stationary  position,  it  drops 
down  rapidly  and  rupture  begins.  In  terms  of  Irwin's  concept  (1958),  the  rate  of  strain 
energy  released  at  this  stage  is  equal  to  a  critical  value,  Gc,  which  marks  a  state  of 
instability. 
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Compression  tests  have  revealed  (Bieniawski  (1967)  that  the  major  type  of 
unstable  fracture  propagation  is  still  following  the  direction  of  major  principal  stress. 

It  has  been  commonly  observed  by  researchers  (Bieniawski,  et.al.  (1969),  Brady  (1970)) 
that  the  unstable  fracturing  process  is  greatly  affected  by  the  type  of  loading,  r-'fe  of 
loading  and  type  of  machine  used  in  the  testing. 

Mechanical  properties  of  rock  containing  flaws  or  cracks  are  somewhat  differ¬ 
ent  from  those  of  intact  rock.  From  a  statistical  point  of  view,  flaws  or  cracks  in  the 
rock  can  be  assumed  to  be  randomly  distributed  and  oriented.  Therefore,  no  definite 
expression  of  the  material  properties  can  be  obtained  unless  the  distribution  function 
of  flaws  is  known.  In  a  series  of  papers  Walsh  (1965-a,b,c)  derived  mathematical 
formulas  for  effective  compressibility,  modulus  of  elasticity  and  Poisson's  ratio  as 
functions  of  mean  crack  length  and  mean  unit  volume  which  contains  the  mean  erack 
length.  Those  mathematical  expressions  were  shown  to  be  qualitatively  justified  by 
a  limited  number  of  experiments.  In  general,  it  was  found  that  compressibility  decreased 
as  the  cracks  closed  under  pressure  and  remained  almost  constant  with  further 

change  in  pressure.  Brace  (1965)  reported  that  linear  compressibility  is  a  function 
of  applied  pressure  and  is  related  quantitatively  to  porosity,  grain  size  and  dimen¬ 
sional  orientation.  If  the  clastic  modulus  is  defined  as  the  tangent  to  the  stress- 
strain  curve,  the  initial  tangent  in  region  I  (Figure  II-4)  represents  Young's  modulus 
of  a  rock  with  initial  open  cracks  (E();  the  constant  slope  in  regions  II  and  III  repre¬ 
sents  the  modulus  for  rock  with  closed  erack  (Ef);  and  the  initial  tangent  to  the  un- 


loading  curve  in  regions  II  and  III  represents  the  intrinsic  modulus  of  an  intact  rock 
(E0)  (Walsh  (1965-b)).  The  magnitudes  of  these  quantities  can  be  ordered  as 

E.  <  E,  <  E 
1  f  o 

Walsh  (1965)  observed  that  the  Poisson's  ratio  of  a  rock  having  open  cracks 
is  slightly  lower  than  that  of  a  solid  rock.  The  value  of  Poisson's  ratio  ranges  from 
zero  to  0.5  as  Ej  changes  from  zero  to  the  value  of  E0t  and  v f  is  greater  than  u0 
but  within  the  limit  of  0.5. 

For  study  of  crack  initiation  and  propagation  in  tensile  stress  fields,  Griffith 
(1920,1924)  proposed  a  theory  postulating  the  propagation  of  elliptical  cracks  as  an 
instability  phenomenon.  A  crack  was  expected  to  extend  if  the  energy  released  in 
propagation  exceeded  the  surface  energy  of  the  additional  surface  created  by  the 
extended  crack.  Initial  flaws  exist  randomly  in  a  solid  body  and  have  random  orienta¬ 
tion.  Fracture  propagation  can  occur  at  the  extremities  of  these  flaws.  The  theory 
is  briefly  outlined  in  Appendix  A  to  this  report.  For  nonuniform  multi-axial  stress 
fields,  a  crack  may  propagate  until  it  stabilizes  on  reaching  a  region  of  low  stress. 
Also  for  finite  systems,  this  effect  of  stress  redistribution  associated  with  crack 
propagation  may  be  significant  in  relation  to  the  extension  of  other  microcracks  in 
the  system. 

Crack  propagation  can  also  occur  at  the  ends  of  closed  cracks  in  compressive 
stress  fields.  McCIintock  and  Walsh  (1962)  developed  the  modified  Griffith  theory 
applicable  to  these  situations.  This  theory  is  also  discussed  in  Appendix  A. 


Extensive  surveys  (Hoek  and  Bieniawski  (1965))  of  published  experimental  data 
have  shown  the  Griffith  and  modified  Griffith  theories  to  be  satisfactory  criteria  for 
fracture  initiation.  Tests  by  Brace  (1963)  and  Lajtai  (1971)  confirmed  that  in  com¬ 
pressive  stress  fields,  cracks  started  at  the  ends  of  closed  joints,  in  the  direction 
normal  to  the  maximum  tensile  stress  and  then  became  parallel  to  the  direction  of 
the  major  compressive  stress  applied  to  the  specimen.  At  this  stage,  the  crack 
stabilized  and  additional  compression  was  needed  for  further  growth.  Brace  (1963) 
observed  cracks  propagating  at  stress  levels  much  lower  than  required  to  extend  a 
single  crack.  Experimental  evidence  (Hahn  (1971))  exists  to  show  that  the  energy 
absorption  associated  with  fracture  may  be  one  or  two  orders  of  magnitude  higher 
than  the  energy  of  surfaces  created.  Moreover,  the  experimental  results  show  con¬ 
siderable  scatter.  Thus  energy  balance  does  not  appear  to  be  a  useful  means  to 
examine  crack  propagation  in  rock. 

Irwin  (1956)  proposed  the  use  of  stress  intensity  factor  to  predict  crack  pro¬ 
pagation.  Stresses  or  displacements  could  be  used  to  obtain  the  value  of  the  stress 
intensity  factor  and  propagation  would  occur  if  the  value  equaled  or  exceeded  a 
critical  value  characteristic  of  the  material.  Compliance  or  energy  balance  methods 
have  also  been  used.  The  finite  element  method  has  been  used  to  evaluate  the  stress 
intensity  factor  for  complex  geometries  and  loadings.  However,  this  approach  is 
unsatisfactory  inasmuch  as  the  extent  of  propagation  and  stress  redistribution  associ 
ated  with  such  propagation  cannot  be  directly  evaluated.  Also,  the  crack  geometry 
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hay  to  ho  known  before  hand.  It  is  not  possible  to  study  initiation  of  fractures  at 
mieroeracks  or  Griffith  Haws  arbitrarily  oriented  and  located  in  the  material. 

in  the  current  research  program,  the  Griffith  and  modified  Griffith  theory  were 
used  to  predict  mioroorack  initiation  and  propagation.  Post-fracture  behavior  across 
a  crack  was  simulated  on  the  assumption  that,  an  open  crack  cannot  transmit  tensile 
ami  shear  stresses,  in  the  finite  element  idealization,  an  element  was  presumed  to 
have  flaws  in  all  directions  such  that  fracture  was  entirely  dependent  upon  the  stress 
field.  Also,  for  small  elements  it  is  reason  Me  to  assume  that  if  an  element  cracks, 
the  fracture  extends  throughout  the  element  at  a  constant  orientation.  With  this 
modeling,  the  exact  location  of  the  'Tack  within  an  element  is  unimportant. 


CHAPTER  m.  THE  FINITE  ELEMENT  METHOD 


3.1.  Basie  Concepts  in  Direct  Method  of  Solution 

A  boundary  value  problem  can  be  stated  u.  the  form 

Au  =  f  on  F  (III-l) 

where  u  is  the  unknown  function  to  be  determined,  A  is  an  operator,  and  f  is  the 
"forcing"  function.  F  is  the  domain  of  interest  and  may  be  an  open,  connected, 
bounded  spatial  region  embedded  in  or  in  a  cartesian  product  x  (.0,  oo  ) 
where  [0,oo)  is  the  non-negative  time  interval.  In  addition  to  the  field  equation 
(III-l) ,  there  will  be  some  conditions  to  be  satisfied  on  boundary  A  of  F.  For  A 
linear  positive,  it  can  be  shown  that  Equation  (III-l)  has  a  unique  solution.  Nec¬ 
essarily,  any  approximate  solution  will  in  general  not  coincide  with  the  unique 
solution  of  Equation  (III-l)  and  consequently  no  approximate  solution  is  expected 
to  satisfy  the  field  equation  as  well  as  the  boundary  conditions  completely. 

Solutions  to  engineering  problems  as  well  as  the  forcing  functions  are,  in 
general, bounded  and  therefore  belong  to  L2 ,  the  space  of  square  integrable 
functions.  L2  is  a  Hilbert  space.  However,  u  may  be  contained  in  a  subset  D 
of  L2  such  that  A  is  defined  on  D.  We  assume  that  D  is  dense  in  L2.  If  the  set 
of  functions  j  0^,  k=  1,2, . . .  ,00  j.  is  a  basis  in  D,  then  any  function  u£  L2  can  be 


expressed  as  an  infinite  sum: 


>■  =  E  \  «k 


k-l 


(in-2) 


A  scheme  to  generate  approximate  solutioi  s  is  to  use  a  finite  set  of  terms  n  the 
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infinite  sum  above.  Thus,  we  use 

N 

“  Z  \  <t>k  (III- 3) 

k  1 

as  an  approximation.  The  approximation  process  then  consists  of  appropriate  choice 
°f  N,  </>k  and  the  coefficients  a^.  Several  alternative  procedures  are  available.  The 
finite  element  method  is  a  special  process  of  selection  of  finite  subset  of  the  basis 
I  1  •  The  coefficients  ak  are  evaluated  by  requiring  the  approximate  solution  to 
satisfy  the  field  equations.  Often  a  more  systematic  approach  is  to  use  a  variational 
formulation  and  obtain  ak  by  requiring  the  approximate  solution  (Eq.  Ill- 3)  to  satisfy 
the  variational  principle.  Ritz'  method,  Galerkin's  method,  least  squares  method, 
all  belong  to  this  category. 

The  finite  element  method  is  well  documented  in  literature  (Zienkiewicz  (1972), 

Bell  and  Holand  (1969),  WPAFB  Conferences  (1965,1968,1971),  Felippa  (1966),  Clough 
(1960,1965)).  Its  theoretical  basis  (Oden  (1969),  de  Arantes  e  Oliveira  (1968),  Zlamal 
(1968),  Melkes  (1970),  Aubin  (1972))  and  relationship  to  variational  principles  (Melosh 
(1963),  Pian  and  Tong  (1969))  have  been  examined.  Essentially  a  finite  element  ideali¬ 
zation  partitions  the  spatial  region  F  into  a  finite  number  of  nontrivial  discrete  elements 

or  subregions.  The  geometry  of  the  elements  is  defined  by  a  set  of  points  in  space  called 
the  nodal  points  of  the  system . 

Over  an  clement  e,  let  an  approximation  to  u  be 

Ne 

=  Z  aek  £ek  (in-4) 

k=l 
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or  in  matrix  form 


e 

u 


{0e}Vi 


(in-5) 


where  {  <pe  }  is  a  row  veetor  consisting  of  0^  as  its  c’cments  and  |  °  |  is  a 

e 

column  vector  of  coefficients  a^  .  Evaluating  the  function  at  nodal  points 


ki 


(III-6) 


where  {  u.  J  is  the  vector  of  nodal  point  values  of  the  function  and  J  is  the  matrix 

T 

of  base  functions  evaluated  at  each  nodal  point.  Rows  and  columns  of  ^  0.  J  are 
linearly  ’ndependent.  If  square,  the  matrix  is  invertible.  If  the  number  of  nodal 
points  is  not  equal  to  the  number  of  independent  base  funetir>ns,  a  least  squares  pro¬ 
cedure  ean  be  used  for  inversion.  Hence,  we  can  write 


lael  -  i[*ef]’1l«?l 

-  [a I1  KI 


(in-7) 


r  -7  c  -i  T 

where  A  -  [<p\  ] 


Substituting  Equation  (III-7)  in  Equation  (III-5) 

T  _  -1 

0°  1  !  A 


e 

u  - 


l  Ui  I 


(in-8) 


i0e}T|uf  1 


where  -  0ej  ean  now  be  regarded  as  a  set  of  interpolating  functions  relating  nodal 


point  values  of  a  function  to  the  value  at  an  arbitrary  point  within  the  element  e. 


3.2.  Variational  Methods  in  Continuum  Mechanics 


3.2.1.  Underlying  Philosophy 
Let  A  be  an  operator  such  that 

A  :  V - -V*  (III— 9) 

where  V,  V*  are  linear  vector  spaces  over  a  field  F  or  a  finite  connected  subset 
thereof.  For  u£V,  we  introduce  the  operator  equation 

Au  =  f  .  3  .  f£V*  (ni-10) 

Let  B  be  a  bilinear  map  on  V* 

B  :  V*  x  V*  S  (in- 11) 

where  S  is  a  linear  vector  space.  To  each  ordered  pair  of  vectors  u,v£V*,  B 
assigns  a  point  B(u,v)£S,  such  that 


B(au1  +  u2,v)  =  aB(u1,v)  +  B(u2,v) 


(in- 12) 

B  (u,  a  Vj  +  v2)  =  a  B(u,  v^)  +  B(u  ,  v2) 

For  simplicity,  in  the  sequel  we  shall  use  the  notation  <  u,v  >  as  equivalent  to  B(u,v). 

In  order  to  set  up  variational  principles  corresponding  to  the  field  problem  ex¬ 
pressed  by  Eq.  (Ill- 10),  we  introduce  a  pair  {ft,  A  } .  The  first  element,  ft  ,  defined 
through  a  bilinear  map  B  on  V*,  is  a  function  of  the  quantities  appearing  in  Eq.  (III-10), 
and  for  given  A  and  f  can  be  regarded  as  a  function  of  u.  The  other  element  of  the  pair, 
A,  is  a  variation  operator  defined  on  the  range  of  ft  in  S.  By  appropriate  selection  of 
B,  ft  ,  A  ,  different  variational  formulations  of  a  problem  are  realized.  In  the  context 
of  obtaining  an  approximate  solution  to  Eq.  (III-10),  if  B  is  continuous  and  V,  S  are 
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metrie  spaces,  the  distance  of  an  approximation  Uq  from  the  correct  solution  u 
is  intuitively  related  to  tho  distance  between  their  respective  images  under  fl  ,  thus 
providing  a  basis  for  study  of  convergence  of  sequences  of  approximate  solutions. 

Many  variational  formulations  use  R,  the  space  of  real  numbers,  as  the  range 
of  B.  For  this  ease,  two  alternative  approaches  are  used.  One  defines  A  so  that 


Af)  >0  (III- 13) 

whenever  (ni-10)  is  satisfied.  This  yields  a  class  of  'minimum'  principles.  Another, 
somewhat  more  versatile  approach  applicable  to  bilinear  maps  in  general  is  to  define 
A  so  that 


A  n  -  0  (III- 14) 

is  equivalent  to  Eq.  (III-10).  An  interesting  feature  of  this  second  alternative  is  that 
if  0 ,  A  are  defined  such  that 

AQ(u)  --  Au-f  (III- 15) 

Eq.  ( III—  1 4)  is  direetly  equivalent  to  Eq.  ( III— 10) .  Thus  A  can  be  viewed  as  a  gradient 
operator  or  derivative.  In  literature,  this  has  >ccn  identified  with  the  Gateaux  deriva¬ 
tive,  the  Frcehet  derivative  or  the  variation  operator  used  by  Gurtin.  The  Gateaux 
derivative  of  n,  denoted  by  G^j  fi(u),  is  defined  such  that 


<  G -  n  (u),  u  > 


(III- 16) 


Here  <  ,  >  indicates  the  bilinear  map  appropriate  to  the  problem  and  u£V*  is  the 
'path'  or  the  'direction'  of  the  derivative  such  that  ufau£V;aisa  scalar.  We  note 
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here  that  for  the  notion  of  a  limit,  we  require  a  suitable  topology  is  S  and  continuity  of 
Q.  1’hese  analytical  aspects  are  not  discussed  in  this  report.  It  is  assumed  that  such 
conditions  where  necessary  arc  satisfied.  'The  right  hand  side  of  Eq.  (Ill- 16), 
denoted  V  17  (u,  !i) ,  is  the  Gateaux  differential  of  17  at  u  in  the  direction  u; 
assumed  to  be  linear  in  u  to  establish  Q(u).  Following  Vainbeig  (1964),  we  write 
a  linear  Gateaux  differential  as  D  17  (u,u).  To  ensure  a  norm,  S  has  traditionally 
been  identified  with  R,  the  set  of  real  numbers,  although  clearly,  this  is  not 
necessary. 

If  V*  and  S  arc  norrned,  the  Freehct  derivative,  denoted  by  F_  17 (u)  is  defined 

u 

such  that 

<  I-'n  17 (u ) ,  u  >|  Jam  I  ft(uiu)  -  17  (u) I  (III-17) 

|u|— ►O 

where  |  |  denotes  a  norm.  Gurtin  (1966,  1964)  regards  17(ui  >iu)  as  a  function  of  a, 

a  scalar.  Then, 

A-  17  (u)  -  -7 -  17  (u  miu) 

u  ua 

a  0  (III- 18) 

provided  the  derivative  exists.  A-  17 (u)  is  expected  to  be  of  the  form 

A_  17 (u)  <  Au  -  f,  u  >  ( III— 19) 

such  that  vanishing  of  the;  variation  for  arbitrary  11C  V*  is  equivalent  to  Eq.  ( ni— 10) .  Indeed 
A_  17(u)  can  be  shown  to  be  identical  to  the  linear  Gateaux  differential.  For  Eq.(III-lO)  to 
be  equivalent  to  Kq.  (Ill- 19), in  conjunction  with  Eq.  (Ill— 14) ,  the  bilinear  map  must  have 
the  property* 

*  We  do  not  use  special  notations  to  distinguish  between  zero  elements  of  V,  V*,  S. 
The  nature  of  the  zero  is  apparent  from  the  context  in  which  it  is  used. 
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<  u,  o  >  (I  <  o,  V  >  ;  v  u,  v£  V* 


(III— 20) 


<  u,  v  >  0  ->  u  o  or  v  0  (IN-21) 

A  class  of  variational  principles  dircetly  uses  Eq.  (Ill— 20)  along  with  its  con¬ 
verse,  i.e.,  if  <  u,  v  >  =  0  V  v£V*,u  -  0.  Thus  if  we  define  fl  as  a  bilinear  map 
on  V*  x  V*  sueh  that 

£l(u,v)  =  <  Au  -  f,  v  >  (HI-22) 

fi(u,  v)  -  0  V  vGV*  direetly  implies  Eq.  (Ill- 10).  This  approach,  for  operators 
with  suitable  symmetry  properties,  leads  to  principal  of  ’virtual  work'  type.  In  this 
formulation,  based  on  the  orthogonality  of  zero  with  the  linear  veetor  spaee  V*,  the 
variation  operator  plays  no  part  and  may  be  taken  to  be  the  identity  operator.  For 
obtaining  approximate  solutions  to  Eq.(IEHO),  Eq.(III-20)  ean  be  used  as  a  starting 
point  for  several  well-known  procedures.  For  example,  if  v  =  Au  -  f,  defining  0 (v)  = 

|  <  v,  v  >|  leads  to  a  generalization  of  the  least  squares  method. 

3.2.2.  A  Formulation  for  Linear  Operators 

Consider,  corresponding  to  the  field  problem  expressed  by  Eq.  (Ill— 10) , 

fi(u)  -  <  u,  Au  >  -  2<u,f  >  (III-23) 

where  <  ,  >  denotes  a  bilinear  map  satisfying  Eqs.  (Ill— 12 ,  III-20,  III-21).  Using 
the  variation  operator  of  Eq.  (Ill— 18) ,  we  obtain  for  linear  A, 

A-  Q(u)  -  <  u,  Au  >  +  <  u,  Au  >  -  2<  u, f  >  (III-24) 

If  A  is  symmetric,  i.e. , 

<u,  Av  >  =  <  v,  Au  >  ;  u,  vGV  (ni-25) 

then  _  _  _ 

O(u)  =  <  Au,  u  >  +  <  u,  Au  >  -  2<  u,f  >  (III-26) 
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If,  further,  the  bilinear  map  is  symmetric /i.  e. , 


<u,v  >  <v,u  >  ;  u,  v£V*  (III— 27) 

it  follows  that 

A-  fl(u)  2<u,  Au  -  f  >  ( III— 2 8) 

Eq.  ( III— 28)  seen  with  Eqs.  (III-20)  and  (III-21)  implies 

A-  Q(u)  -  0  for  arbitrary  u£V*  (III-29) 

if  and  only  if  Eq.  (Ill— 10)  is  satisfied. 

We  recognize  A-  17  (u)  as  the  generalized  Gateaux  differential  such  that 
Eq.  (Ill— 2 8)  implies 


Gj  Q(u)  -  2  (Au  -  f)  ( III— 30) 

Clearly  G-  Q(u)  vanishes  if  and  only  if  Eq.  ( III— 10)  is  satisfied. 

The  foregoing  discussion  shows  that  for  linear  bounded  symmetric  operator  A, 
solution  of  equation  Au  -  f  is  equivalent  to  vanishing  of  variation  of  ft( u)  where  in  the 
pair  {  Q,  A  }  ,  fl  =  <"  u,  Au  >  -  2  <  u,f  >  ,  and  A  is  the  variation  operator  defined 
by  Eq.  (III-18)  or,  alternatively,  is  the  Gateaux  derivative.  The  bilinear  map  <  ,  > 
satisfied  Eqs.  (III-12), (III-20) , (ill-2 1)  and  (III-271. 

If  the  symmetric  bilinear  map  has  its  range  in  R,  and  in  addition  to  Eqs.  <111-12), 
(III— 20) ,  and  (III— 27)  satisfies 

<u,  u>  >  0,  V  u  i  0  (III— 31) 

the  bilinear  map  defines  an  inner  product.  In  this  ease,  for  A  positive,  i.e. 


<  u,  Au  >  >  0,  u  /  0  (III— 32) 


#Magri  (1972)  has  shown  that  for  every  linear  operator,  A,  there  exists  a  bi¬ 
linear  form  B  such  that  A  is  symmetric  in  the  sense  of  B. 
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using  the  pair  {  P. ,  A}  where  fi(u)  is  defined  by  Eq.  (Ill— 23)  and  A  is  defined  as 

A-  fi(u)  =  ft(u  +  u)  -  Q(u)  (III- 3 3) 

minimum  principles  can  be  developed  for  positive  or  positive  definite  operators.  A 
comprehensive  exposition  has  been  presented  by  Mikhlin  (1905)  for  this  case. 

Traditionally,  variational  principles  have  been  based  on  an  inner  produet, 
which  is  a  bilinear  functional  satisfying  Eq.  (III-31)  but  not  Eq.  ( III— 2 1) .  The  motiva¬ 
tion  has  been  to  use  positive  property  (Eq.  (III-32))of  the  operator  to  establish  unique¬ 
ness  of  the  solution.  Also  the  functionals  have  their  range  in  R,  the  set  of  reals,  a 
Hilbert  space.  However,  other  bilinear  maps  are  available. 

The  above  discussion  for  a  single  field  variable  u  is  easily  extended  to  the 
ease  of  several  dependent  field  variables.  If  there  are  n  variables,  V  is  defined 
as  the  direct  sum  space  Vj  0  V2  0  ...  0  Vn  and  an  element  u£V  is  an  n-tuple 
(UpUg,  . . . ,  un)  with  Uj  6  Vj  for  i  =  1,2, . . .  ,n.  A  bilinear  map  on  V  is  now  defined 
as 

<  u, v  >  =  <  Ui, v.  >  i  <  u2,v2  >  !■  . . .  •  <  un , vn  >  (III- 34) 

12  n 

(no  sum  on  n) 

where  <  ,  >  is  a  bilinear  map  from  Vr*  x  V*  »  S.  Symmetry  of  each  of 
<,  >r  implies  symmetry  of  <u,v>.  The  operator  A  is  defined  as  a  two-dimensional 
array  of  linear  operators  with  n  x  n  elements  such  that  an  element  A-  maps  a  set 
M.j  C  Vj  into  V.  .  Corresponding  to  Eq.  (III-25),  we  define  symmetry  of  A  as 

<  uj  ,  Ajj  Uj  >  =  <  Uj,  Ajj  U|>  (no  sums)  (III-35) 
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For  diagonal  elements  of  the  array,  Eq.  (III-35)  is  identical  to  Eq.  (III-25).  The 


operator  equation  (III- 10)  is  now  a  set  of  linear  equations  for  the  problem 


Ajj  Uj  -  ft  ,  i,j  =  1,2,...  ,n 


(III- 36) 

Equation  (III-23)  will  then  define  ft(u)  and  we  shall  require  each  of  the  variations 


A  u.  ft  (u)  =  2  <  uit  Ay  Uj  -  f{  >  (III- 37) 

for  ujEVi;  i  =  1,2, . . .  ,n  to  vanish.  Examples  of  this  type  of  formulation  have  been 
presented  by  Sandhu  and  Pister  (1970, 197  1).  Often  problems  which  are  not  in  linear 
symmetric  form  can  be  manipulated  to  write  them  in  the  form  of  Eq.  (Ill— 36)  with 
Ay  satisfying  Eq.  (Ill- 35).  Symmetry  of  the  operator  matrix  leads  to  extended  vari¬ 
ational  principles  such  that  the  unique  intersection  of  the  sets  of  solutions  associated 
with  these  alternative  formulations  is  the  problem  solution.  Generalizations  to  include 
nonlinear  operators  on  the  diagonal  of  the  operator  matrix  arise  as  natural  extensions. 

These  aspects  of  the  problem  have  been  discussed  in  detail  by  Sandhu  and  Pister  (197 ’,) 
who  also  presented  ageneral  discussion  of  variational  principles  in  continuum  mechanics  (1972). 


3.2.3.  The  Elastostatics  Problem 

The  variational  formulation  of  problems  in  elasticity  has  been  discussed  by 
Washizu  (1968)  and  by  Sandhu  and  Pister  (1971).  Herein  we  only  present  generalized 
potential  energy  formulation  in  which  the  variation  of  the  functional 

£1  =  /  <2U  -  Uj  u1JfJ  -  2  <ij  "-Jj  ♦  Uj  j  Tjj  -  2  V,  >  dF 


\f  “i  <*ij  "j  - 2  ‘t>  ds  -  ,  /<“ i  -  2V  <rij  "j  d' 
X>\  A 2 


(III— 38) 
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vanishes  at  the  correct  solution  to  the  field  equations 


(in-39) 


-  U 

U  (1,3) 


( III— 40) 


(in-4i) 


subject  to  the  boundary  conditions 


'■«  nj  =  'i  on  A 

u  =  u  on 
i  i  2 


(III— 42) 


(III- 4  3) 


Here^lt^2  are  complementary  subsets  of  j$,the  boundary  of  F.  U  is  the  strain  energy 

function,  Ojj  ,  are  respectively  the  components  of  the  symmetric  Cauchy  stress  tensor 
and  the  infinitesimal  strain  tensor.  Uj,  arc  the  components  of  the  displacement  vector 

and  the  body  force  vector,  nj  are  direction  cosines  of  the  outward  normal  to^S 

Using  the  symmetry  property  of  the  field  equations,  and  requiring  uj  to  identi¬ 
cally  satisfy  the  strain  displacement  equation  as  well  as  the  boundary  conditions  on 
^2*  ^cb  (ni-38)  reduces  to  the  potential  energy  functional 


^2  ~  f  (U  -  uj  f{)  d  F  -  /  Uj  tj  d  s 


(III— 44) 


To  allow  for  initial  stresses.the  Eq.  (III-40)  is  written  as 


dU 

& i-  =  - - -  t  (T. . 

«  d'ij  1J 


(III— 45) 


Equation  (III-46)  is  the  functional  used  in  development  of  stiffness  type  finite 
element  formulations. 

In  the  analysis  for  progressive  failure  according  to  Griffith's  theory,  the 
above  formulation  was  used  assuming  the  system  to  be  piecewise  linear  elastic  — 
changes  in  structural  properties  being  introduced  in  a  stepwise  fashion  as  elements 
in  the  system  develop  cracks  in  a  sequential  manner. 

3.2.4.  Variational  Formulation  for  Elastic- Plastic  Solids 

The  field  equations  of  incremental  plasticity  in  symmetric  form  are 


(m-47) 

Here  E^i  are  components  of  the  symmetric  tensor  relating  the  components  of  the 

incremental  stress  tensor  cr..  to  the  components  of  the  incremental  strain  tensor  i 

J  ij 

“i  »  Fk »  °"ij  are  respectively  the  components  of  the  incremental  displacement  vector, 
the  body  force  vector  (including  seepage  forces  if  any),  and  the  initial  stress  tensor 
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for  the  specific  increment  under  consideration.  Following  Sandhu  and  Pister  (1970, 
197 1)  the  governing  functional  is 


J2  =  J  (U  -  Uj  Fi  +  Uijj  o^)  dF  -  J  Ui  tt  d  s 
F 


(m-52 


which  is  of  the  same  form  as  fig  in  Eq.  (HI-46)  except  for  the  displacement  components 
uj  being  replaced  by  the  incremental  displacement  uj  and  U  replacing  U.  Explicitly, 

o 

introducing  the  seepage  force  into  the  formulation, 


=  p  L  +  n  . 


(III-53 


where  n  is  the  hydrostatic  pore  water  pressure.  Eq.  (III-52)  then  becomes 


J3  /(U-Ujpft-Ui  ni+<xi] 
F 

f  •  - 

J  Uj  t|  ds 


«Tii)  t,F 


(III— 54 


The  above  formulation  would  reduce  to  that  given  by  Washizu  (1968)  if  the  equilibrium 
equation  were  to  be  written  in  the  form 


(<r ..  .  +  <r ..  .)  +  (F.  +  F.  )  -  0 
ij,)  iJ.J  i  i 


(III— 55 


with  the  boundary  condition  ti  -  ti  +  tj 


(III-56 


Then  upon  using  symmetry  and  the  specialization  of  displacements  satisfying  the 
boundary  conditions  and  the  strain  displacement  equation,  the  governing  functional 


would  have  the  form 


J4  " 


/(U-i  F)  dF  -  /  u,  tj 


(III— 57 


The  finite  el  rnient  procedure  for  elastic-plastic  Mohr-Coulomb  rock  described  in 


Chapter  IV  was  based  on  the  variational  principle  using  the  functional  J„. 


3.3.  Matrix  Formulation  of  Field  Equations 

Following  the  procedure  indicated  in  seetion  3. 1,  a  finite  dimensional  approxi¬ 
mation  to  the  problem  solution  is  written  using  the  values  of  the  function  at  the  nodal 
points  of  the  system  as  the  generalized  coordinates.  The  spatial  domain  F  is  parti¬ 
tioned  into  a  finite  number  of  nontrivial  disjoint  elements  sueh  that  their  union 
approximates  F.  Eaeh  point  in  F  belongs  to  one  distinet  element.  The  ambiguity 
regarding  points  on  interelement  boundaries  is  resolved  by  ensuring  that  the  function 
value  will  be  the  same  regardless  of  whieh  element  the  point  is  assigned  to.  This 
is  the  requirement  of  'compatibility'  on  finite  element  modelling.  For  a  point  within 
an  element,  the  function  value  is  expressed  in  terms  of  nodal  point  values  (the  gen¬ 
eralized  coordinates)  through  a  set  of  interpolating  functions.  The  approximate 
function  is  inserted  in  the  governing  functional  so  that  the  functional  can  now  be 
regarded  as  a  function  of  the  generalized  coordinates.  Vanishing  of  variation  of 
the  functional  yields  the  matrix  equations  for  evaluation  of  the  generalized  coordinates 
and henee  the  approximate  solution.  We  give  below  the  formulation  for  ineremental 
analysis  of  elastie-plastie  solids.  The  treatment  for  the  elastic  ease  is  similar. 

The  functional  Jg  in  Equation  (III— 54)  has  increments  of  the  displacement  vector 
as  the  field  variable.  In  an  approximate  solution,  the  ineremental  displacement  field 


within  the  mth  element  is  defined  in  terms  of  the  nodal  point  values  as 

•  m  ,  ,m  i 1  * 

u.  (x)  {0  1  u. ; 


(III— 5  8) 


■  m 


Here  u.  (x  )  is  the  ith  component  of  the  incremental  displacement  vector  at  spatial 


loeation  defined  by  the  veetor  x  in  a  global  reference  frame,  |0m}  is  the  set  of 
displacement  interpolation  iunetions  and  -1  u  j }  is  the  set  of  the  ith  components  of 
nodal  point  incremental  displacements  for  the  entire  system. 

The  strain  displacement  relationship  ean  be  written  as 


(III— 59 ) 


Here  -jem  (x)  j  is  the  redueed  strain  tensor.  For  the  two-dimensional  ease, 


jem<x)) 


rxv 


T 


(in-60) 


r  ml 

L^e  J  *s  transformation  matrix  derived  from  the  displacement  interpolation 
functions  -[0mj  by  suitable  differentiation  and  rearrangement  of  terms  and  (uj 
is  the  veetor  of  nodal  point  incremental  displacements  for  the  system. 

Writing  the  stress-strain  relationship  as  matrix  [Hm]  for  the  nith  element  and 
replacing  integration  over  the  spatial  domain  by  sum  of  integrals  over  individual  ele¬ 
ments,  the  governing  functional  Eq.  (111-54)  becomes 


where 


[K]  = 
l«l)  = 


iu|T[K]|ul+2.u)T{M1}  -2 


E 

m-l 


(III— 61) 
(III— 62 ) 

(III- 6 3) 
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(in- 64) 


2, 

m  =  1 


f[<pm]  (pH  dvm 


E  fir]  [rjT  >rj  dSm 

ml 


(III- 65) 


Application  of  the  variational  principle  to  Equation  (III-61)  yields  the  matrix  equation 

[k]  {  u }  =  {r|  (III— 66) 

where 

{r1  =  {M1|  -  ^M2)  +  |  P  j  (ffl-67) 

{  M  |  ,  -  p}  respectively  represent  the  contribution  of  the  residual  stresses, 

the  body  forces  and  the  boundary  loads.  For  incremental  formulation, the  residual 
stresses  are  the  stresses  at  the  beginning  of  the  increment.  Pore  water  pressures 
are  included  in  the  form  of  seepage  force  in  the  body  force  contribution. 

For  the  elastic  ease,  the  function  to  be  approximated  is  the  displacement  vector. 

•  • 

Equations  (III— 66) ,  (III— 67)  will  apply  with  the  difference  that  i  u  ' ,  -  R  ’  will  be  re¬ 
placed  by  u  |  ,  <  R  and  j  q-j  j  will  represent  the  initial  stresses  in  the  system. 

The  computer  programs  developed  both  for  elastie-plastic  rock  and  progressive 
failures  following  Griffith  theory  used  a  quadrilateral  element  made  from  four  constant 
strain  triangles.  This  element  as  well  as  other  quadrilateral  elements  have  been  dis¬ 
cussed  by  Wilson  (1965)  and  Dougherty  et.al.  (1968). 


5i 
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3.4.  Initial  Stresses  in  Rock 


The  initial  stresses  present  in  rock  in  the  'insitu'  state  arc  the  stresses  corre¬ 
sponding  to  reference  displacement  field  generally  taken  as  zero.  These  are  easily 
included  in  the  stress-strain  relationships  by  assuming  that  the  total  stress  is  the 
sum  of  the  initial  stresses  and  the  stresses  associated  with  the  displacements.  Thus 

"ij  =  Eijkl  fkl+*ij  <m-G8) 

where  ojj ,  f ^  are  components  of  stress  and  strain  increments  and  o-^  are  components 
of  the  initial  stress.  With  this  formulation,  the  initial  stresses  appear  in  the  variational 
formulation  and  vector  |  in  Equation  (III-61)  represents  a  pseudo-load  corresponding 
to  the  initial  stresses.  This  technique  is  the  basis  of  stress-relief  methods  used  by 
various  investigators.  The  analyses  described  in  Chapters  IV  and  V  use  this  approach. 

3.5.  Incremental  Procedures 

In  finite  clement  analysis,  allowance  for  incremental  loading  or  incremental  con¬ 
struction  is  quite  straightforward.  For  any  step  of  loading  or  construction  (excavation), 
the  initial  state  is  used  as  a  reference  state  and  increments  of  stresses  and  displace¬ 
ments  worked  out  for  the  particular  step.  Thus 

Wn  Mn  ■  Mn  I111’69' 

where  [.K]n,  j  Au|n,  jAR^,  are  respectively  the  system  stiffness  matrix,  the 
incremental  displacement  vector  and  the  incremental  loading  vector  for  the  nth  step 
of  construction /excavation/loading  In  eases  where  [k]r  is  dependent  upon  displace¬ 
ments  or  stresses,  iterative  procedures  or  Runge-Kutta  methods  can  be  used  to  im- 


prove  accuracy.  In  the  research  described  in  this  report,  essentially  the  method 
was  ’initial  stiffness'  for  each  increment  in  the  case  of  elastic-plastic  analysis. 
For  the  case  of  progressive  cracking,  the  influence  of  stress  redistrioution  and 
stiffness  change  due  to  cracking  was  allowed  for  in  an  iterative  process. 
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CHAPTER  IV.  FINITE  ELEMENT  PLANE  STRAIN 
ANALYSIS  OF  ELASTIC- PLASTIC  MOHR-COULOMB  ROCK 

4.1.  Review  of  Previous  Work 

Attempts  to  use  the  finite  element  method  for  stress  find  deformation  analysis 
of  elastie-plastie  materials  may  be  classified  on  the  basis  of  the  modeling  of  mechan¬ 
ical  behavior  or  the  computational  technique.  Early  attempts  noticed  the  nonlinearity 
of  deformation  response  to  applied  loads.  The  material  was  regarded  as  nonlinear 
elastic  having  stress  or  strain  dependent  elastie  'moduli'.  In  order  to  extend  the 
results  of  the  me-dimensional  test  to  the  six-dimensional  stress-state,  an  equivalent 
stress  to  equivalent  strain  curve  was  adopted.  Bilinear  (Wilson  1963)  or  multilinear 
(Zienkiewicz  1967)  approximations,  Richard- Goldberg  law  (1965),  Ramberg-Osgood 
law  (1943),  Koudner's  hyperbolic  equation  (1963),  cubic  spline  "unctions  (Desai  1971) 
have  been  used  to  approximate  test  data  by  smooth  or  piecewise  smooth  curves  for 
ease  of  data  handling  within  the  computer. 

Swedlow  and  Yang  (1965)  used  the  normality  rule  and  Drueker's  (1952)  method 
for  evaluating  the  constant  of  proportionality  in  the  Prandtl-Reuss  equation.  This 
was  a  tangent  modulus  approach.  It  did  not  satisfy  the  condition  f  =0. 

M areal  (1968)  and  Mareal  and  King  (1967)  used  similar  formulation.  Reyes  (1966) 
and  Reyes  and  Deere  (1966)  used  a  rate  of  work  equation  to  develop  stress-strain 
relations  for  Mohr-Coulomb  materials  under  plane  strain  conditions.  This  was 
consistent  with  the  incremental  theory  of  plasticity.  Felippa  (1966)  developed 
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equations  for  incremental  theory  of  plasticity  which  include  Reyes  work  as  a  speciali¬ 
zation.  Marcal  and  King  (1967)  and  Zienkiewicz  et  al.  (1969)  use  a  formulation  iden¬ 
tical  to  Felippa's.  Baker  et  al.  (1969)  used  Reyes  formulation  in  their  work.  Yamada 
(1968)  used  a  rate  of  work  in  distortion  to  set  up  equations  for  von  Mises  materials. 
Isakson  et  al.(19f7)  developed  the  incremental  stress-strain  relationship  for  kinema¬ 
tic  hardening  using  Ziegler's  (1959)  modification  of  Prager's  rule  (1955). 

Numerical  procedures  for  elastic-plastic  analysis  have  been  either  the  incre¬ 
mental  type  using  Euler  or  Runge-Kutta  methods,  or  the  initial  tangent  type  using  the 
initial  stress  or  initial  strain  approach.  In  the  first  type  the  system  stiffness  matrix 
has  to  be  developed  at  each  increment  whereas  in  the  initial  tangent  method,  the  same 
stiffness  is  used  throughout  with  the  effect  of  nonlinearity  being  introduced  as  a  cor¬ 
rective  pseudo-load  in  conjunction  with  an  iterative  solution  scheme.  The  initial  strain 
approach  was  used  by  Argyris  (1965),  Gallagher  (1362),  Lansing  et  al.  (1965)  among 
others.  Zienkiewicz  (1969)  used  the  initial  stress  approach.  Baker  et  al.  (1969) 
found  the  iterative  procedure  to  be  unsatisfactory  as  convergence  was  often  very  slow. 

All  the  research  workers  mentioned  above  used  triangular  elements  in  their 
analysis.  It  is  well-known  that  the  triangular  elements  do  not  give  correct  stress 
field.  In  this  report  quadrilateral  elements  were  used  to  get  a  good  stress  distribution. 
In  the  present  research  program  the  incremental  approach  reflected  in  the  variational 
formulation  of  the  problem  was  used.  The  formulation  follows  Felippa  (1966)  and 


Reyes  (1966). 


nwarunmujuiuiPiiwiTgpwa 


4.2.  jtr ess-Strain  Relations  for  Mohr-Coulomb  Materials 

The  Mohr-Coulomb  failure  criterion  for  isotropic  materials  is 

Tf  "  c  +  °n  tan0  (IV— 1) 

where  Tf ,  <rn  are  the  failure  shearing  stress  and  the  normal  stress  on  the  failure  plane 
and  c,  <p  are,  respectively,  the  cohesion  and  the  angle  of  in  ernal  friction  for  the  mate¬ 
rial.  Drucker  and  Prager  (1952)  proposed  a  generalization  of  the  Mohr-Coulomb  law 
to  a  cone,  in  the  three-dimensional  principal  stress-space,  symmetrical  about  the 
principal  diagonal.  The  loading  surface  is  defined  as 

j 

f  =  q  J f  t-  J22  -  k  -  0  (IV-2) 

where  J^,  are  the  first  and  the  second  invariants  of  the  stress  tensor,  and  a,k  arc 
material  constants. 

Using  Drucker  and  Prager's  generalization  of  the  Mohr-Coulomb  law,  Reyes 
(1966)  developed  the  stress-strain  relationship  for  elastic  perfectly-plastic  solids 
under  plane  strain  conditions.  A  rate  of  work  equality  was  used.  The  same  relation¬ 
ship  can  be  obtained  b j  directly  using  Equation  (IV-2)  in  Felippa's  general  approach 
already  discussed  in  Chapter  II. 

4.2.1.  Stress- Strain  Formulation  for  Elastic- Plastic  Mohr-Coulomb  Solids 
Define  a  tensor  with  components  q..  such  that 

“D  *  T7T  <IV-3> 

For  Drucker  and  Prager's  generalization  of  Mohr-Coulomb  yield  surface 


~  T*-*s*y**ifc.  -----  ^ 


f  =  a  Jj  +  J22  =  k.  This  gives 


q..  =  a  R..  +  - *- 

11  11  j 

V 


af„  t- 

IJ 


—  (cr..  -  I  J  f...) 

r  A  li  3  1  lJ 


(IV— 4) 


where  s.j  are  components  of  the  stress  deviation  tensor  and  8^  is  Kronecker's  delta. 


Define  a  tensor  with  components  p..  such  that 

6  f 

Pij  =  - 7T 

3  *'ij 


(IV— 5) 


For  perfect  plasticity,  the  yield  surface  is  independent  of  the  plastic  strain  and  p^  =  0. 
Let  q  be  a  reduced  representation  for  q„  such  that 


q  '  <  qll  q22  q33  2ql2  Zq23  2q31  > 


(IV-6) 


where 


qll  '  a  +  -7T  <*11  -  \  Jl> 


q22  =  a  +  (”-22  "  A  Jl> 


a  +  — 1 —  (or  -  1  J  ) 

.  1  '33  3  1 

J  2  2 


r  .1  12 

d22 


(IV— 7) 


Following  the  notation  used  in  Chapter  II,  Equation  (11-30^, 
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H  =  1  +  9  a2  A 
G 

vt  K,  G  are  respectively,  the  Poisson's  ratio,  the  bulk  modulus  and  the  shear  modulus 
in  elasticity. 

For  plane  strain 

•  • 

a  -  hQ  e  +  h.  tr  a  e  „ 

2  ‘aa  1  ap  eap 

b  =  h  i  +  h  tr  f 
1  *aa  3  a{3  *a(3 

where  summation  on  repeated  indices  is  over  the  range  1,2. 

4. 3.  Analysis  of  Progressive  Fai'ure 
4.3.1.  Basic  Methodology 

In  applying  the  finite  element  method  to  elastic-plastic  continuum,  it  was 
assumed  that  the  stress  field  is  constant  within  each  element.  Under  applied  loads, 
an  element  was  assumed  to  be  either  totally  yielded  or  elastic.  For  sufficiently 
small  elements,  this  appears  to  be  a  reasonable  assumption.  The  stiffness  of  the 
system  represents  the  collective  stiffness  of  its  elements.  For  each  element  after 
yield,  the  stress-strain  behavior  was  defined  by  the  relationships  developed  in  section 
4.2.  Thus  as  each  element  yields,  the  system  stiffness  is  altered  resulting  in  a  non¬ 
linear  structural  response.  In  the  study  of  progressive  failure,  it  is  important  to 
allow  for  the  effect  of  this  nonlinearity.  In  the  current  research  program,  the  system 
was  assumed  to  be  stepwise  linear  within  yield  of  successive  elements.  In  the  incre- 


mental  procedure,  the  load  was  divided  into  increments  with  each  increment  large 
enough  to  introduce  yield  in  another  element.  This  was  subject  to  checks  on  the 
validity  of  assumption  of  linear  stress-strain  relationship.  As  several  elements 
may  yield  at  almost  the  same  load,  a  lower  limit  on  load  increment  size  was 
used  to  avoid  excessive  computational  effort.  Also,  elements  yielding  within  a 
prescribed  load  range  of  the  computed  increment  for  any  step  were  treated  as  if 
these  had  yielded  simultaneously  with  the  elements  controlling  the  increment  size. 

The  total  load  vector  j  R  |  was  treated  as  a  sum  of  load  increments  j  ARj  j> ,  i  =  1, 

2, . . . ,  n  where  n  is  the  total  number  of  increments,  not  necessarily  equal. 

{  R}=  £  {  ARi  }  (W-Il) 

The  first  increment  of  load,  j  AR^|  applied  to  an  elastic  system  was  suffi¬ 
ciently  large  to  allow  at  least  one  element  to  reach  the  yield  surface.  To  evaluate 
the  load  increment  for  any  step  as  a  proportion  of  the  load  yet  to  be  applied  to  the 
system,  a  stress  ratio  was  introduced.  Figure  IV- 1  shows  the  calculation  of  the 
stress  ratio  Sr  for  a  typical  element.  The  point  A  represents  the  initial  stress  state 
and  C  the~5tress  state  that  would  result  if  all  the  load  was  applied.  The  curve  f  -  k  =  0 
represents  the  yield  surface.  If  A  and  C  both  were  in  the  interior  of  the  yield  surface, 
all  the  load  could  be  applied  and  Sr  was  defined  as  unity.  Howe -er,  if  A  was  in  the 
interior  and  C  in  the  extarior  of  the  yield  surface,  clearly  it  would  not  be  possible  to 


Ww,.  . . . 
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load  the  element  to  C  and  that  it  could  only  be  loaded  to  B  before  onset  of  yield.  Here  B 
is  the  intersection  of  line  AC  with  the  yield  surface.  Then 


I A  B| 

I A  C  | 


(IV-12) 


Calculation  of  Sr  is  simple, 
with  A  and  C  respectively. 


Let  (cr^).,  (°ij)f  represent  the  stress  states  associated 
Then  the  stress  state  for  B  is  o-..  such  that 

i] 


aSj  =  <Vi  (1'Sr>  +Fr<*il>f  <IV-13> 

Sr  is  calculated  from  the  relationship 

f  (°i j )  -k  =  0  (IV- 14) 

The  stress  ratio  Sr  for  an  element  represented  the  fraction  by  which  the  load  would 
have  to  be  scaled  to  ensure  that  after  application  of  the  load  increment,  the  element 
was  within  or  on  the  yield  surface.  The  value  of  Sr  was  calculated  for  all  the  elements. 
The  element  with  the  lowest  stress  ratio  was  the  next  to  yield  and  governed  the  system 
load  ratio.  Sj.  thus  defined  the  load  increments  in  progressive  failure  assuming  step¬ 
wise  linear  behavior.  As  each  load  increment  was  applied,  cumulative  stresses  and 
displacements  were  calculated.  The  stresses  at  the  end  of  the  ith  step  were  the  initial 
stresses  for  the  (i  +  l)th  step. 

4.3.2.  Nonmonotonic  Loading 

To  allow  for  non-monotonic  loading,  at  any  structural  increment  or  excavation, 
a  pilot  analysis  is  carried  out  assuming  all  elements  to  be  linearly  elastic.  In  case  of 
unloading,  any  element  already  on  yield  surface  can  unload  either  plastically  or  elastic- 


ally  as  shown  in  Figure  IV- 2  .  The  two  modes  are  distinguished  by  the  fact  that  in 
elastic  loading,  the  strain  decreases  whereas  if  the  element  stays  on  the  yield  sur¬ 
face  in  the  unstable  region  the  strain  will  increase.  If  the  element  is  assigned  the 
stress-strain  behavior  corresponding  to  the  yield  state,  it  is  forced  to  stay  in  the 
yield  surface  because  the  stress-strain  laws  satisfy  f  -=  0.  In  the  pilot  analysis, 
assuming  elastic  behavior, the  element  will  cither  follow  the  path  AB  or  AC.  Path 
AC  corresponds  to  elastic  unloading.  Path  AB  would  indicate  increasing  strain  under 
applied  load  and  therefore  the  element  is  likely  to  stay  in  the  yield  surface.  This  pre¬ 
liminary  information  is  used  to  assign  appropriate  material  behavior  to  each  element 
previously  in  state  of  plastic  yield. 

4.3.3.  Allowance  for  Nonlinear  Stress-Strain  Behavior 
The  stress-strain  relationship  for  plasticity  is  stress-dependent.  As  the 
stresses  change  during  application  of  a  load  increment,  it  is  desirable  to  reflect 
this  change  in  the  stress  computation.  Assuming  cjj  cjj  ^  to  be  the  compliances 
at  the  initial  and  final  stress  states,  respectively,  we  may  write 

Ut<0  =  11,(0)  +  iiL-LIU.®  *  f, 

2  J 

where  ,  u^0)  ,  a  Pj  are  the  final  displacement,  the  initial  displacement  and  the 
increment  load  vectors,  respectively.  Further 

cij(f)  =  cij  (°kl(f)  ) 


and 


°kl(f)  =  °-kl(ui(f)> 
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STRAIN,  £ 


FIG.  IV-2.  Loading  or  Unloading  for  an  Initial  State  of  Plastic  Yield 
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An  iterative  scheme  would  consist  of  the  following  steps 
ut(f)("+1>  =  Ui(°)  +  cij(0)  +  cij(f)(n)  A 


<rki(f)(n)  =  (rkl  (ui(1)(n)) 


n  =  0,  1, 


Such  an  iterative  procedure  has  been  used  by  Sandhu  (197  3)  along  with  a  criterion  for 
convergence  which  defines  the  maximum  increment  of  a  Pj  consistent  with  the  non¬ 
linear  characteristics  of  cij^.  However,  repeated  solution  for  displacements  may 
be  very  expensive  in  terms  of  computation.  Also,  because,  in  general,  tne  plasticity 
will  be  confined  to  local  regions,  the  displacement  solution  is  not  likely  to  change 
significantly.  Therefore,  it  is  considered  sufficient  to  assume  constant  displacements 


and  to  iterate  only  on  stress  using  mean  stiffness  k^j  as  under: 
o-JfMn+l)  =  ^.(0)  +  kiJkI(0)  +kijkl(f)(n) 


A  €kl 


kijkl  _  ^ijkl  (^mn  ) 

where  a  is  the  strain  increment.  The  procedure  converges  very  rapidly.  The 

procedure  is  shown  schematically  in  Figure  IV- 3a. 

An  alternative  procedure  is  to  use 

o-  Mn+1>  =  «r  <°)  +  k  (m)(n)  Ae 

ij  ^ij  +  ijkl  ^€kl 


where  the  mean  stiffness  corresponding  to  the  load  increment  is  defined  as  the  stiff¬ 
ness  at  the  mean  stress  for  the  increment,  i.e. 


65 


Ae  =  STRAIN  INCREMENT 


FIG.  IV-Ha.  Calculation  of  Stress  lor  a  Strain  Increment  -  Method  1 


k  (m)(n) 
Kijkl 


Cijkl  ( 


(O' 

°mn  +  o- 


mn 


(f)(n) 


This  approach  is  illustrated  in  Figure  IV-3b.  Both  the  alternatives  were  found  to  be 
satisfactory. 

4.4.  Examples  of  Application 

A  computer  program  for  plane  strain  analysis  of  elastic-plastic  Mohr-Coulomb 

rocks  was  developed  using  the  mathematical  model  and  the  discretization  procedures 
described  in  the  foregoing  sections.  The  technique  was  used  to  solve  several  problems. 

There  are  very  few  problems  in  the  theory  of  plasticity  for  which  closed  form  solutions 

are  available.  However  Naghdi  (1957)  solved  the  problem  of  an  clastic  perfectly  plastic 

wedge  under  uniform  loading  on  one  face  (Fig.  IV-4).  The  Naghdi  solution  is  for  a 

wedge  infinite  in  extent,  made  of  von  Mises  material  and  loaded  in  plane  strain.  This 

type  of  material  is  a  special  case  of  Mohr-Coulomb  material  obtained  by  setting  the 

angle  of  internal  friction  0=0.  This  example  was  used  earlier  by  Baker  et  al.  (1969) 

to  verify  their  computer  code.  Another  example  Involves  analysis  of  stresses  and 

deformation  of  a  notched  bar  of  perfectly  plastic  von  Mises  material.  The  last  example 

is  of  frequently  used  laboratory  test.  Excellent  agreement  with  theoretical  results  for 

the  wedge  were  obtained. 

4.4,1.  Elastic-Plastic  Wedge 

Figure  (IV-4)  shows  a  finite  element  idealization  for  the  Naghdi  wedge.  Figures 
(IV-5)  and  (IV-6;  show  the  theoretical  and  computed  results  for  the  distribution  of 
radial  and  circumferential  stresses  at  various  stages  of  loadings.  The  angle  V  de¬ 
notes  the  angle  upto  which  the  yielding  has  progressed  from  the  boundaries.  Figure 
(IV -7)  shows  the  radial  strain  distribution  at  various  ranges. 
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Naghdi's  theoretical 
solution 


A  Finite  element  results 


Generally  the  agreement  between  results  computed  by  the  method  outlined  and 
the  exact  analysis  was  found  to  be  good. 

4.4.2.  Notched  Specimen 

A  notched  specimen  of  perfectly  plastic  von  Mises  material  with  a  90°  notch 
and  subjected  to  a  uniformly  distributed  load  was  analyzed.  The  finite  element  ideali¬ 
zation  for  one  quarter  of  the  specimen  is  shown  in  Figure  IV-8.  The  value  of  c  was 
taken  as  12. 15  kilogiam  per  square  millimeter  (corresponding  to  a  yield  stress  of 
24.3  kg.  /mm.  in  uniaxial  tension  test)  and  0  was  set  equal  to  zero.  Figure  IV-9 
shows  the  principal  stresses  at  a  load  intensity  of  19.2  kilogram  per  square  milli¬ 
meter.  Figure  IV-10  shows  contours  of  failure  ratio  for  this  load  as  well  as  the  boun¬ 
dary  of  the  plastic  zone  at  different  values  of  the  load.  This  problem  was  solved  by 
Mareal  and  King  (1968)  and  Zierkiewiez  et  al.  (1969).  The  results  are  slightly  differ¬ 
ent  from  those  presented  in  this  report  for  reasons  discussed  earlier. 


FIG.  IV-8.  NOTCHED  SPECIMEN—  FINITE  ELEMENT  IDEALIZATION 
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SCALE:  1  —  100  kg/mm2 


FIG,  IV-9.  Notched  Specimen  (One  Quarter)  -  Elastic- Plastic  Stress  Distribution  Under 
Uniform  Load  of  19.2  kg. /mm,^  Applied  at  Ends 
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Dashed  line  show 
extent  of  region  «n 
plastic  yield 


fV-10.  Notched  Specimen  -  Contours  of  Failure  Ratio  at  Uniform 
Extensional  Load  of  ID. 2  kg. /mm.  Applied  at  Ends 
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CHAPTER  V.  FINITE  ELEMENT  PLANE  STRAIN  ANALYSES  OF 
PROGRESSIVE  CRACKING  OF  ROCK  FOLLOWING  GRIFFITH'S  THEORY 


5.1.  Review  of  Previous  Work 

Stability  of  excavations  in  rock  is  deeply  influenced  by  the  presence  of  cracks 
and  fissures.  These  discontinuities  may  be  preexisting  or  might  arise  as  a  conse¬ 
quence  of  the  stress-redistribution  associated  with  excavation.  Classical  methods 
of  analysis  are  inadequate  for  study  of  initiation  and  propagation  of  fracture. 

Several  attempts  have  been  made  to  apply  the  finite  element  method  to  jointed 
rock  systems.  Zienkiewicz  et  cl.  (1968)  proposed  the  'no  tension'  analysis  pro¬ 
cedure.  This  consisted  of  the  following  steps: 

1.  Analysis  of  the  system  treating  intact  rock  as 
linear  elastic,  isotropic. 

2.  Check  to  identify  tensile  principal  stresses,  if 
any,  in  various  elements. 

3.  Reanalysis  assigning  zero  resistance  to  defor¬ 
mation  in  the  direction  of  the  principal  tensile 
stresses. 

4.  Repeat  2  and  3  to  convergence,  i.e,  until  the 
solution  shows  no  appreciable  tension  anywhere. 

In  order  to  economize  on  computational  effort,  a  stress  relief  procedure  was  intro¬ 
duced.  The  tensile  principal  stresses  were  relieved  by  introducing  equivalent  nodal 
point  loads  using  the  linear  elastic  stiffness  L'or  an  iterative  correction  scheme.  The 
procedure  has  poor  convergence  characteristics  and  even  after  several  cycles  of  itera¬ 
tions,  tension  zones  do  not  completely  disappear.  Moreover,  for  the  case  of  one 
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principal  stress  being  tensile  and  the  other  compressive,  relieving  tensile  stress 


by  applying  equivalent  nodal  point  loads  amounts  to  using  a  non-symmetrical  consti¬ 
tutive  relationship.  In  the  present  research  effort,  a  correction  was  introduced 
in  the  no-tension  approach  to  correctly  simulate  the  orthotropic  material  behavior 
of  cracked  elements.  For  intact  linear  isotropic  elastic  elements, 


C  C 
11  12 

C21  C22 


(V-l) 


where  o-j,  <r2  are  the  principal  stresses;  <2  are  the  principal  strains  and 

CH»  ci2*  C2i*  C22  are  elements  in  3  symmetric  matrix.  In  the  no-tension  method 
^1  >  °»  °2  <  °»  °i  would  be  relieved  by  applying  nodal  loads £bTj|<^1j  whereJbTj 
represents  the  equilibrium  transformation.  The  stress  i  is  retained  i.e,  for  the 

Li 

cracked  element, 


C21  C22 


(V-2) 


is  the  stress- strain  relationship.  Clearly,  for  symmetry  in  material  behavior,  the 
part  of  <r2  corresponding  to  C2l  must  also  be  relieved  so  that  for  cracked  element 


(V-3) 


Even  the  corrected  approach  was  not  considered  worthwhile.  Aside  from  the  poor  con¬ 
vergence  of  the  numerical  solution  procedures,  there  is  a  basic  objection  to  eliminating 


tension  from  several  elements  simultaneously.  Cracking  must  necessarily  be  progres- 
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sive  and  as  each  eraek  forms,  the  rock  in  the  immediate  vicinity  experiences  a  signifi 
cant  stress-redistribution.  The  no-tension  approach  does  not  allow  an  ..',o  sequential 
cracking  of  elements  and  is  therefore  not  realistic. 

Duncan  an-!  Goodman  (1068)  examined  the  effect  of  preexisting  joints  at  arbitrary 
orientations.  A  linear  elastic  analysis  was  used  to  establish  the  stress  field.  Normal 
and  shear  etre  <^s  on  arbitrarily  oriented  planes  (ubiquitous  joint)  were  computed  and 
the  shearing  strength  based  on  Mohr-Coulonib  law  compared  with  the  shearing  stress. 
This  analysis  is  useful  in  predicting  local  failure.  However,  the  analysis  did  not  take 
into  account  the  stress-redistribution  and  progressive  failure  associated  with  local 
failure. 

In  a  comprehensive  report,  Duncan  and  Goodman  (1968)  also  used  an  orthotropie 
continuum  to  simulate  rock  with  orthogonal  sets  of  parallel  and  evenly  spaced  joints. 
Again  the  progressive  failure  arid  deformation  could  not  be  allowed  for. 

For  preexisting  joints,  it  is  possible  to  simulate  their  mechanical  behavior 
through  the  use  of  fwo-dimensional  elements  or  one-dimensional  elements.  The  two- 
dimensional  elements  have  the  drawback  of  poor  'aspect  ratios'  (Duncan  and  Goodman 
(1968))  in  the  ease  of  very  thin  joints  leading  to  inaccuracy  in  the  results.  One- 
dimensional  elements  to  simulate  joints  were  developed  by  Goodman  et  al.  (1968) 
and  have  nroved  quite  useful.  These  elements  can  transmit  shear  as  well  as  com¬ 
pressive  stresses.  Figure  ( V—  1)  shows  the  representation  discussed  by  Dunean  and 


Goodman. 


Propagation  of  preexisting  cracks  has  also  been  studied  in  conjunction  with 
the  concept  of  stress-intensity  factors.  Assuming  the  crack  geometry  to  be  known, 
finite  element  procedures  were  developed  by  Ch;m  et  al.  ( 1 1)7 0)  and  Gross  et  al. 
(19f>8)  to  determine  the  stress  intensity  factors  at  crack  tip.  Displacements,  stresses 
and  compliance  were  all  used  as  the  basis  for  calculation  of  the  factors,  l'o  improve 
accuracy,  Wilson  (1971),  Byskov  (1970),  Levy  (1971)  introduced  stress-singularity 
elements  it  eraektips  (Figs.  V-2  &  V-3).  Piaxi(197  1)  suggested  use  of  hybrid  elements. 
The  stress-singularity  elements  improve  the  accuracy.  However,  they  are  severely 
restrictive  in  the  study  of  crack  propagation  as  the  analysis  applies  only  for  crack 
tip  at  the  cente-  "f  the  elemcn  .  Use  of  quadrilateral  elements  with  stress  criterion 
is  sufficient,  in  most  cases,  to  establish  the  stress-intensity  factors. 

A  common  drawback  of  all  the  previous  wurk  is  that  the  geometry  of  the  cracks 
and  their  location  must  be  known  beforehand.  For  new  cracks  originating  at  Griffith 
flaws,  the  crack  geometry  is  not  known  beforehand.  The  present  research  program 
treated  the  problem  of  fracture  initiation  at  arbitrarily  oriented  Griffith  flaw's  using 
the  stress  formulation.  Assigning  orthotropic  no-tension  material  behavior  to  cracked 
elements,  it  was  poss’blc  to  trace  the  progressive  failure  of  rock  following  Griffith's 
theory.  The  procedural  details  arc  given  in  the  following  sections. 

5.2.  Analysis  for  Crack  Init  ition 

In  the  ease  of  preexisting  cracks,  it  is  sometimes  possible  to  plar  a  finite  ele¬ 
ment  mesh  including  the  crack  surface  as  a  free  boundary.  However,  wl  ere  intact 
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rock  cracks  under  changes  in  stress  environment,  it  is  necessary  to  cheek  for  crack 
initiation  in  accordance  vith  Griffith  or  modified  Griffith  theory.  As  element  pro¬ 
perties  change  only  at  cracking,  for  a  given  initial  mechanical  state,  a  linear  clastic 
analysis  is  valid  upto  cracking  of  an  additional  element  in  the  system. 

Assuming  Griffith  flaws  are  oriented  in  every  direction  <n  each  element,  an 
analysis  based  on  stresses  in  the  element  is  carried  out  to  verify  crack  initiation. 

The  theoretical  basis  for  this  is  outlined  in  Appendix  A.  Upon  application  of  full  load, 
several  elements  may  satisfy  the  criterion  for  crack  initiation.  In  this  case,  a 
stress  ratio  is  established  for  each  element.  The  stroSS  ratio  is  the  factor  by 
which  the  stress  increment  must  be  multiplied  such  that  the  total  stress  corresponds 
to  the  critical  stress  environment  for  the  element.  The  stress  ratio  is  a  function 
of  the  initial  stress  and  the  stress  path.  The  minimum  stresB  ratio  indicates  the 
next  element  to  crack  in  a  progressive  failure  sequence. 

5.3.  Analysis  of  Progressive  Cracking  of  Rock 

As  each  element  eraeks,  its  stiffness  ehanges  reflecting  a  change  in  the  system 
stiffness.  Thus  the  total  load  oeformation  behavior  ».f  rock  is  piecewise  linear  in  a 
finite  element  representation,  changes  in  slope  being  associated  with  sequential  crack¬ 
ing  of  various  elements.  In  actual  situations,  the  nonline,  rity  may  be  continuous. 

In  the  process  of  analysis,  the  system  is  assumed  to  be  piecewise  linear  elastic. 
At  each  stage  the  live  load  is  applied  and  the  str  *ss  response  calculated.  As  fracture 
of  one  element,  next  in  the  sequenec,  marks  a  ehangc  in  stiffness,  stress-ratio  is 
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ealeulated  for  eaeh  element  as  shown  in  Figure  (V-4).  Point  A  represents  the  initial 
state  of  stress.  Point  C  along  stress  path  AC  represents  the  stress  state  if  the  system 
were  linear  elastie  for  the  load  increment  applied.  If  point  C  satisfied  the  fracture 
criterion,  it  was  neeessary  to  find  the  stress  ratio  Sr=  sueh  that  the  point  B 

defined  a  eritieal  stress  state.  As  fraeture  angle  depends  upon  the  stress  environ¬ 
ment,  calculation  of  the  stress  ratio  involved  an  iterative  procedure.  For  fraeture 
angle  0  corresponding  to  C,  a  stress  ratio  was  calculated.  Then  for  this  stress  ratio  , 
the  corresponding  stress  state  was  evaluated  and  the  eritieal  angle  for  this  stress  state 
represented  a  better  approximation  to  the  eorreet  fraeture  angle.  Convergence  was  r  ipid. 
The  stress  ratios  for  all  elements  were  compared.  The  minimum  value  represented 
the  end  of  the  particular  step  in  the  multilinear  stress  path  dependent  system,  and  the 
element  yielding  the  minimum  ratio  was  the  nextto  eraek.  Repetitive  application  of 
this  procedure  defined  the  progressive  failure  of  rock. 

In  order  to  ensure  stability  in  computation,  the  total  load  was  applied  in  several 
increments.  This  served  to  keep  the  stress  ehange  in  individual  elements  for  a  given 
load  increment  within  reasonable  limits.  Figure  (V-5)  illustrates  incremental  loading 
analysis.  For  a  typieal  load  increment  AP,  the  initial  stiffness  yields  point  C  and 
allowance  for  sequential  fraeture  of  several  elements  within  this  load  increment  shifts 
the  (P, u)  point  to  B.  The  different  arrows  indicate  sequential  eorreetions  to  disp»  ce¬ 
ment  allowing  for  progressive  fraeture. 
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Note:  p  calculated  for  state  B  will  in  general  not  be  the 
same  as  Q  for  state  C. 


FIG.  /-4.  Calculation  of  Stress  Ratio,  S 


-  Initial  Stress  for  the  Load  Increment 
C  =  Stress  for  Full  Load  Increment 
B  -  Scaled  Increase  cf  Stress  to  Crack  Initiation 


5.4.  Modelling  of  Cracked  Rock 


Referred  to  axes  of  material  symmetry,  the  constitutive  relationship  for  linear 


orthotropic  isothermal  elasticity  can  be  written  as 


'V 

1 

eT 

~v2 

E2 

~v3 

E3 

0 

0 

0 

' 

'll 

*22 

~V1 

E1 

1 

*2 

~v3 

E3 

0 

0 

0 

<r 

22 

*33 

'  = 

^1 

E1 

E2 

1 

0 

0 

0 

< 

'33  > 

y23 

0 

0 

0 

G23 

0 

0 

*23 

y31 

0 

0 

0 

0 

G31 

0 

a3l 

.  yi2  , 

0 

0 

0 

0 

0 

<*12^ 

K  '12  . 

Symmetry  of  the  compliance  matrix  requires 

E1  E2  E3 


(V->) 


( V— 2) 


For  plane  strain  *33  =  0,  ^23=  ^31“  ® 


( V—  3) 


Using  Eq.  (V-i), 

'r33  =  *'3<1  <V-4> 

Substituting  Equation  (V-4)  in  Equation  (V-l), 
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Inverting 


(1-  i*j  c3)  (1-  1^3)  -  k  ^22  ( 1  +  ^ 3) 2 


El<1"^2^)  El‘,2<1U,3) 


(V-G) 


tL-, 

where  k  — L 

e2 


Upon  development  of  fracture  surface  in  an  element,  the  plane  of  the  crack  can 


be  regarded  as  a  plane  of  mechanical  symmetry.  Also  the  crack  plane  is  a  principal 


plane  for  the  element.  This  involves  a».  „ssumption  that  the  crack  is  planar  and  ex¬ 


tends  throughout  the  element.  It  appears  to  be  reasonable  for  sufficiently  small 


element  size. 


Assuming  that  cr^  corresponds  to  the  stress  normal  to  the  fracture  plane,  the 


elastic  modulus  is  reduced  to  a  very  small  value.  The  element  is  still  capable  of 


withstanding  stresses  <r  22  parallel  to  the  crack  plane.  However,  when  <r2 2  also  attains 


values  such  that  another  fracture  plane  develops  within  an  element,  the  material  can¬ 


not  sustain  any  load.  This  is  modelled  by  letting  E2  also  become  very  small. 


As  crack  orientation  may  not  be  along  the  reference  directions,  it  is  necessary 


to  transform  th,_  stress-strain  relationship  from  the  principal  axes  to  the  reference 
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axes.  Let  the  tingle  of  rotation  be  (3  .  Then  from  Figure  (A -4) 

(3  y  -  a 

where  Y  the  angle  between  c rx  and  o-j 

a  the  angle  between  cr2  and  the  fracture  plane 

For  Griffith  tailure  criterion,  it  can  be  shown  that 

a  -  tan 


1  |  (k-1)  sin20  *  1  -  y  (k2  -  1)  sii\2d  »  1  j 

A 


\  (1  -  k)  sin  26  j  -0 

where  k  -  cr2  /  o-j  ,  and  Q  is  identical  to  (3  defined  by  A-4  in  Appendix  A. 

For  the  modified  Griffith  criterion, 

n  f) 

Q  =  T  -6 

where  6  <s  the  same  angle  as  (3  defined  by  Eq.  (A-15)  in  Appendix  A. 
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L*t  the  relation  between  principal  strains  and  strains  in  global  coordinates  be 
represented  by 
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Here  j  j  ,  |  |  are  the  strains  referred  to  the  principal  axes  and  the  reference 

axes  respectively  and 
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The  relation  between  principal  stress  and  stresses  in  the  reference  frame  is 


It  can  be  shown  that  if 
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5.5.  Pre-existing  Discontinuities 

Pre-existing  discontinuities  in  rock  may  be  initial  weak  planes  or  initial  open 
joints.  In  either  case,  strength  of  material  would  be  reduced.  For  an  initial  open 
joint,  the  orthotropic  stress-strain  relation  described  in  Section  5.4.  can  be  used,  and 
the  initial  fracture  angle  would  be  defined  by  the  initial  crack  plane  or  by  the  inclination 
of  joint.  In  the  case  of  initial  weak  plane,  we  propose  an  orthotropic  stress-strain  rela¬ 
tion  with  certain  amount  of  shear  resistance  along  the  weak  plane. 

A  relation  similar  to  those  in  Section  5.4  can  be  written  as 
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Following  a  procedure  similar  to  that  in  Section  5.4,  we  obtain 


where  J' 


L  J 

M 


h  •  h- 

r*  i  - 


D' 


<1 

K 


-  <  -  sin  2(3  sin  2(3  cos  2(3  > 


In  general,  the  fracture  plane  does  not  coincide  with  the  weak  plane.  The  elasti¬ 
city  matrix  in  such  a  ease  can  be  obtained  by  deleting  the  third  column  and  the  third  row 
of  [D’3  •  and  in  the  transformation,  3  is  replaced  by  the  angle  between  the  fracture 
plane  and  the  weak  plane. 

To  deal  with  closure  of  open  joints,  it  was  assumed  that  the  crack  opening  is  planar 
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and  that  it  can  be  measured  normal  to  the  crack  plane.  This  measured  quantity  can 
be  either  the  initial  opening  of  a  preexisting  joint  or  the  equivalent  opening  tolerance 
for  an  element  fractured  under  a  load  increment  to  close  again.  The  initial  opening 
is  equivalent  to  a  set  of  nodal  displacements  corresponding  to  the  initial  opening. 

I  he  equivalent  nodal  displacements  thus  obtained  were  further  tr:msformcd  into  cle¬ 
ment  strains,  and  the  possibility  of  closing  was  checked  on  the  basis  of  computed 
strains. 

5 .  G .  Incremental  Excavation 

Execution  ot  an  excavation  project  is  a  sequential  process.  In  discretized 
analysis  procedure,  the  results  at  the  end  of  one  stage  constitute  the  initial  state 
for  the  next  step.  History  of  the  system  in  terms  of  stresses  and  deformations  is 
determined  corresponding  to  the  discrete  steps  in  excavation  or  construction. 

For  each  stage,  an  incremental  loading  analysis  was  used.  The  equilibrium 
equation  corresponding  to  each  load  increment  can  be  written  as 


where  i  denotes  the  ith  increment.  To  take  account  of  the  effect  of  crack  propagation, 
[Kj]  was  modified  for  elements  cracked  and  fApi|  corrected  to  include  the  load 
resulting  from  the  releasing  of  the  initial  stresses  from  the  cracks.  The  total  displace¬ 
ment  and  loading  for  N  increments  at  jth  stage  arc  given  by 
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The  final  displacement  and  loading  for  a  complete  v.  'rstruetion  of  M  stages  are 
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j  =  l  i  1 


5.7.  Examples  of 


The  procedure  developed,  for  analysis  of  progressive  failure  of  roek  following 
Griffith's  or  modified  Griffith  theory,  in  accordance  with  the  mathematical  model 
described  in  the  preceding  paragraphs  were  used  to  solve  several  problems.  On 
cracking  of  concrete  beams,  considerable  volume  of  data  is  available.  The  pro¬ 
cedures  were  also  applied  to  analysis  of  eraek  propagation  in  tunnels  of  circular 
and  elliptical  shape.  Effect  of  incremental  construetion/excavation  was  allowed 
for.  A  simple  model  was  used  to  simulate  loading  in  underground  blasting. 

5.7.1.  Crack  Propagation  in  Concrete  Beams 

Application  of  the  finite  element  method  in  the  study  of  cracking  in  reinforced 
concrete  members  was  first  proposed  by  Ngo  and  Scordelis  in  1967.  In  their  analysis, 
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two-dimensional  triangular  elements  were  used  for  idealization  of  concrete  and  steel 
reinlorcement.  1  he  bond  aetion  between  concrete  and  steel  bars  was  modeled  by 
special  links.  Geometry  of  crack  was  defined  by  disconnecting  nodes  as  preexisting 
openings  and  progressive  craeking  in  concrete  was  not  allowed.  Nilson  (1908)  incor¬ 
porated  an  incremental  loading  process  along  '*ath  a  simple  craeking  criterion  allow¬ 
ing  for  progressive  craeking.  Again,  the  eraeks  were  defined  by  disconnecting  nodes 
when  the  average  stress  at  that  nodal  point  satisfied  the  craeking  criterion. 

rIsing  the  method  of  analysis  described  in  preceding  sections,  sequential 
craeking  in  simply  supported  beams  was  investigated.  Results  of  these  investiga¬ 
tions  arc  summarized  in  the  following  paragraphs. 

A.  A  Plain  Concrete  Beam 

(V-G)  gives  the  configuration  of  a  simply  supported  beam  without  any 
reinforcement.  Crack  commenced  at  midspan  and  propagated  throughout  the  eross- 
seetion.  A  collapse  mechanism  developed  when  the  eraek  extended  to  the  top  element. 
Severe  compression  caused  by  the  hinge  aetion  produces  two  secondary  eraeks  which 
follow  the  modified  Griffith  criterion. 

B.  A  Concrete  Beam  with  Tension  Reinforcement 

A  finite  element  representation  of  a  reinforeed  concrete  beam  with  19  feet 
effective  span  and  a  depth  of  32  inches  is  shown  in  Fig.  (V-7a).  Two  eoneentrated 
loads  of  300  lbs.  each  were  applied  at  two  one-third  points  along  the  total  span. 

The  sequence  of  eraeking  is  indicated  by  numbers  as  shown  in  Fig.  (V-7b)  and  (V-7e). 
The  craeking  may  be  roughly  grouped  into  three  stages.  Upon  application  of  the  load, 
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primary  cracks  start  and  propagate  in  the  middle  third  of  the  span  and  stop 
after  travelling  eertain  distances.  Secondary  eraeks  then  follow  and  spread 
over  the  outer  one-third  of  the  span.  Some  bond  slippage  associated  with  extension 
of  the  primary  eraeks  occurs  at  the  last  stage.  The  load-defleetion  eurve  in  Fig. 

(V-7e)  reflects  the  continuous  proeess  of  cracking. 

C.  A  Concrete  Beam  with  Tension  Reinforcement  (Bresler  and  Scordelis,  1963) 

A  series  of  tests  on  reinforced  eonerete  beams  was  conducted  by  Bresler  and 

Seordelis  (1963).  The  present  method  of  analysis  was  applied  to  two  of  these  beams 
to  predict  the  cracking  behavior. 

Fig.  (V-K)  shows  the  configuration  of  beam  OA-2  and  A-2  in  the  aforementioned 
tests.  All  the  dimensions  and  material  properties  used  are  taken  from  the  paper  by 
Bresler  and  Seordelis  (1963)  (See  Table  V-l,  V-2). 

The  cracking  sequence  for  load  up  to  45  kips  is  shown  in  Fig.  (V— 9) ,  and  the 
load-midspan  deflection  eurve  is  given  in  Fig.  (V-lO).  Agreement  with  the  experimental 
data  both  in  the  eraek  patterns  and  the  load-defleetion  curves  is  excellent.  Referring  to 
Fig.  ( V— 10) ,  increasing  deviation  of  the  eurve  obtained  from  the  finite  element  analysis 
from  the  curve  given  by  tests  would  be  expected  with  further  increase  in  load. 

D.  A  Concrete  Beam  with  Tension  .Compression  and  Web  Reinforcement 

In  addition  to  the  tension  reinforcement  in  beam  OA-2,  compression  as  v  ’.i  as 
web  reinforcement  was  used  in  beam  A-2.  Again  the  cracking  sequence  and  eraek  pat¬ 
tern  were  in  agreement  with  experimental  data.  The  midspan  defleetion  eurve  obtained 
at  a  total  load  of  62  kips  was  slightly  low  r  than  the  eurve  given  by  the  test,  but  in 
general,  agreement  was  excellent  (see  Figs.  V-ll,  V-12). 
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FIG.  V-S.  Finite  Element  Idealization  of  Reinforced  Concrete  Beams  OA-2  and  A-2 


TABLE  v-2.  PROPERTIES  OF  STEEL  REINFORCING  BARS 


AREA  sq. 
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5.7.2.  Progressive  Cracking  Around  a  Semi-Circular  Tunnel 


The  configuration  of  a  lined  semi-circular  tunnel  analyzed  by  Ztenkiewicz 
et  al.  (1968)  using  ihe  no-tension  approach  is  shown  in  Fig.  (V-13).  Limitation 
of  this  approach  has  been  discussed  in  Section  5.1.  rig.  (V-15)  shows  the  initial 
tension  zones  in  the  linear  elastic  solution.  This  is  jimilar  to  the  solution  obtained 
by  Ztenkiewicz  (Fig.  (V-14a)).  Upon  application  of  fracture  criteria,  the  crack  was 
found  to  initiate  and  propagate  sequentially  (Fig.  7-16).  Figs.  (V-16,  V-'7,  V-18) 
show  that  stress  redistribution  around  the  tunnel  resulted  in  reduction  of  the  tensile 
zones.  The  cracks  oecur  mainly  along  the  contaet  between  the  lining  and  the  roek. 
Cracking  of  this  kind  may  be  interpreted  as  failure  of  bond  between  the  concrete 
lining  and  the  roek. 

5.7.3.  Progressive  Cracking  Around  an  Elliptic  Tunnel 

One  major  concern  during  an  incremental  excavation  in  rock  mass  is  frag¬ 
mentation  of  the  material  caused  by  progressive  fracture  occurring  in  the  vicinity 
of  excavated  area.  To  p’  event  such  situation  from  occurring,  supports  or  lining 
ean  be  applied  immediately  following  the  excavation.  A  system  consisting  of  an 
elliptic  tunnel  was  used  to  study  the  progressive  cracking  in  an  incremental  excava¬ 
tion  proeess. 

A.  An  Unlined  Tunnel 

Fig.  (V-19a)  shows  a  system  assuming  excavation  to  be  completed  in  a  single 
step.  Fracture  first  commences  in  the  surfaee  of  the  tunnel,  then  proceeds  upward 
layer  by  layer.  The  cracking  sequence  is  shown  in  ^ig.  (V-19b). 
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Zienkiewics's  Solution  of  Tunnel  Problem 


FIG.  V-16  Sequence  ot  Cracking  and  Final  Tensile  Zones  in  Rock 
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Stress  Scale.  PS i.  *  *192* 

(Arrow  for  Tension) 


FIG.  V-18  Final  Stress  Distribution  Around  thn  Scmi-cin'ular  Tunnel 
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a.  System  Analyzed  b.  Sequence  of  Cracking 

FIG.  V-19  Progressive  Cracking  Around  an  Eliptic  Tunnel 
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It  is  noticed  that  the  cracking  sequence  and  the  patterns  of  crack  are  not  unique 
for  a  system.  In  fact,  they  arc  greatly  influenced  by  the  material  properties.  For 
example,  if<rt,  the  tensile  strength,  is  the  only  varying  parameter  in  the  analysis, 
material  having  a  higher  value  of  would  have  less  extensive  cracking.  Similarly, 
if  p,  the  internal  coefficient  of  friction,  is  changed,  the  cracking  sequence  and  the 
eraek  patterns  will  be  different. 

B.  A  Lined  Tunnel 

The  system  described  in  the  preceding  section  was  reanalyzed  in  three  steps 
using  p  0.8  and  trt  25  psi  (see  Fig,  (V-20)).  Step  one  simulated  excavation  in  a 
single  step  as  in  case  A.  Craeks  appeared  in  the  side  wall  of  the  tunnel  as  predicted 
by  modified  Griffith  criterion  and  ft  the  bottom  by  Griffith  criterion.  This  was  be¬ 
cause  of  severe  compressive  stressts  in  the  side  wall  while  high  tensile  stresses 
occur  at  the  tunnel  invert  and  erown.  Hatches  in  Fig.  (V-20b)  show  the  orientation 
of  the  cracks.  In  step  two,  concrete  linirg  was  introduced.  This  provided  the  tunnel 
wall  with  some  support  due  to  the  tendency  of  the  concrete  lining  to  defleet  outwards 
under  its  own  weight.  The  support  effect  is  increased  if  increase  in  concrete 
temperature  is  allowed  for  or  if  the  concrete/rock  contact  is  pressure  grouted.  A 
region  near  the  horizontal  diameter  developed  double  cracks  in  the  first  step. 

In  the  second  step,  one  set  of  eraeks  in  this  region  closed.  No  further  cracking 
occurred  at  the  end  of  this  step  (see  Fig.  !V-20e)).  In  the  third  step,  a  pressure  of 
2,000  psf  was  applied  on  the  ground  surface  (Fig.  (V-20d)).  New  eraeks  appeared 
at  the  base  of  the  concrete  lining,  propagating  downward  at  an  angle  of  about  45  degrees 
to  the  right.  Meanwhile  the  double  erack  appeared  again  near  the  horizontal  diameter. 
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c.  Crack  Orientations  After  d.  Crack  Orientations  After 

Placement  of  Lining  Application  of  Surface  Loading 

FIG.  V- 20  Incremental  Analyses  of  an  Elliptic  Tunnel 


5-7-4-  Progressive  Cracking  in  Blasting 

The  mechanism  of  breakage  or  fragmentation  of  a  real  rock  caused  by  blasting 
is  very  complicated,  yet  the  major  mode  of  breakage  is  by  brittle  fracture.  When  a 
charge  detonates  in  a  borehole,  blasting  wave  will  propagate  outward  with  very  high 
velocity  and  pressure.  The  magnitude  .and  direction  of  the  velocity  and  pressure  de¬ 
pend  greatly  on  the  properties  of  the  charge  and  the  roek.  They  also  depend  on  the 
detonation  process  and  the  borehole  geometry.  Calculation  of  these  quantities  is  by 
no  means  simple  and  various  theories  as  well  as  empirical  formulas  have  been 
proposed  (e. g.  Brown,  1956). 

Due  to  the  radial  overflow  of  material  accompanying  the  blastinH  wave,  the 
pressure  in  both  tangential  raid  radial  directions  will  decrease  while  radial  cracks 
appear  and  propagate.  Interaction  between  radial  cracks  and  the  tensile  stress  wave 
reflected  by  free  surface  may  increase  the  tensile  stress  at  the  tip  of  those  cracks 
which  are  parallel  to  the  curved  wave  front.  Tests  on  plexiglass  mentioned  by 
Pcrsson  et.al.  (1970)  concluded  that  cracks  propagating  in  a  direction  at  an  angle 
40  to  80  degrees  to  the  normal  of  the  free  surfaee  hare  a  greater  propagation  velocity. 

A  system  shown  In  Fig.  (V-21)  was  analyzed  as  to  simulate  a  bench  type  blast¬ 
ing  .  Inertial  effects  were  not  considered  in  this  analysis,  nor  was  decay  of  pressure 
due  to  gas  entering  the  cracks  allowed  for.  The  ratio  of  the  depth  of  rock  to  the  dla- 

meter  of  the  hole  was  20.  The  effect  of  detonation  of  a  charge  was  simulated  as  sud¬ 
den  application  of  a  radial  pressure  on  the  perimeter  of  the  hole.  Three  different 

values  of  the  pressure  intensity  viz.,  2000,  9000,  4000  psi,  were  considered.  The 
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Finite  Element  Idealization  of 
Circular  Tunnel 
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cracking  sequences  obtained  are  shown  in  Figs.  (V-22),  (V-23)  and  (V-24).  At  2000 


psi,  only  small  radial  eraeks  form  near  the  top  and  the  bottom  of  the  hole.  For  pres¬ 
sure  intensity  equal  to  3000  psi,  cracking  was  severe  around  the  hole,  and  three  major 
cracks  extended  outward.  Finally,  for  pressure  equal  to  4000  psi,  the  pressure  front 
was  pushed  out  further  and  radial  fragmentation  increased  drastically  with  two  major 
eraeks  propagating  to  the  boundary.  Inclination  of  these  two  major  cracks  to  the  nor¬ 
mal  of  the  free  surface  was  30  degrees  and  70  degrees  respectively. 


FIG.  V-22.  Sequential  Cracking  Due  to  Internal 

Pressure  on  Tunnel  Surface  -  2000  psi 
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FIG.  V-24  Sequential  Cracking  Due  to  Internal 

Pressure  on  Tunnel  Surface  -  4000  psi 
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CMAPTKR  VI.  A  PARAMETRIC  STUDY  OF  STRESSES 
IN  STEEL  SU  P1H)RTS  FOR  A  TUNNEL 


(i .  1 .  Statement  of  Work 

The*  comp'itor  program  for  analysis  of  stresses  and  deformations  in  nonhomo- 
geneous  roek  assuming  elastie-plastie  behavior  was  used  in  a  parametric*  study  of 
stresses  in  steel  supports  for  a  tunnel.  Data  for  the  problem  were  provided  by  the 
United  States  Bureau  of  Mines.  Figures  VI- 1  to  VI-4  show  the  configuration  of  the 
tunnel  opening  and  the  four  different  blocking  systems  studied.  Figure  VI— 5  shows 
the  steel  support  structure.  The  initial  stress  field  was  specified  as  hydrostatic 
pressure  corresponding  to  an  overburden  depth  of  1,000  feet  and  a  material  density 
of  165  pounds  per  cubic  foot.  The  objective  of  the  parametric  study  was  to  determine 
the  influence  of  Young's  modulus,  Poisson’s  ratio,  cohesion  and  angle  of  internal 
friction,  upon  the  moments  and  stresses  in  the  steel  supports.  The  range  of  para¬ 
meter  values  specified  by  the  sponsor  is  shown  in  Table  VI- 1. 

6.2.  Method  of  Solution 

6.2.1.  Mechanism  of  Load  Development 

When  a  tunnel  is  excavated,  the  load  carried  by  the  material  removed  must  lie 
carried  by  the  roek  in  the  tunnel  walls  and  by  the  unexeavated  rock  ahead  of  the  face. 
Continued  excavation  at  the  face  results  in  'loss  of  support'  further  increasing  the  stre 
in  the  walls.  For  linear  homogeneous  isotropic  clastic  roek,  it  has  been  shown  (Abel, 
1967)  that  this  effect  is  felt  only  in  a  region  one  diameter  away  from  the  face.  If  sup¬ 
ports  arc  installed  immediately  after  cxeasation,  they  will  share  in  this  transfer  of 
load  as  the  face  is  advanced.  Theoretically,  after  the  excavation  has  progressed  one 
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PARAMETERS 

MATERIAL 

(a)  Rock 

(b)  Shotcrete 

(e)  Steel 

(d)  Timber 

Young’s  Modulus, 

E  (psii 

1  x  10(i  -  10  x  10(i 

2.71  x  10f) 

30  x  10f) 

1.5  x  10(! 

Poisson’s  Ratio, 

V 

0.  1  -  0.5 

0.1 

0.231 

0.03 

Cohesion,  e  (psi) 

1000  -  5000 

3300 

— 

— 

Tan  <p 

(<p  Anule  of  Internal 
Friction) 

1.0  -  2.  1 

1.0 

— 

— 

Material  Density 
(pef) 

105. 0 

150.0 

400.0 

27.0 

(a)  Data  furnished  by  sponsor. 


(b)  "Manual  of  Concrete  Practice,"  Part  2,  19G8,  American  Concrete 
Institute 

E  =  wl.5  33VfT" 
where  w  150  lb. /c.  ft. 

f'c  compressive  strength,  assumed  as  2000  psi 

(c)  "Manual  of  Steel  Construction,"  AISC 

(d)  "Wood  Handbook,"  Forest  Products  Laboratory,  Forest  Research. 
The  properties  arc  for  western  white  pine. 


Table  VI- 1.  Material  Properties 


diameter  ahead  of  the  support,  there  would  be  no  further  development  of  load  on  sup¬ 
ports.  Actually,  it  is  observed  (Abel,  1967)  that  rock  n  ovement  continues  for  a  long 
time  before  reaching  stabilization.  This  results  in  continued  growth  in  the  load  trans¬ 
ferred  to  the  supports.  Also,  exposu)  o  atmosphere,  loss  of  gouge  material,  and 
blasting  damage  may  alter  the  mechanical  properties  of  the  rock  mass  resulting  in 
increasing  deformation  and  increasing  support  stresses.  In  summary,  the  load  develop¬ 
ment  may  be  associated  with  one  or  more  of  the  following  mechanisms: 

a.  Upon  continued  excavation  at  the  face,  the  removal  of  rock  results  in 
increased  rock  load  being  transferred  to  the  walls  of  the  portion  already  excavated 
and  the  supports  in  that  portion. 

b.  Time-dependent  deformation  cf  rock  is  resisted  by  the  supports  resulting 
in  their  taking  on  increasing  load. 

c.  Change  in  material  properties  after  installation  of  supports  will  result  in 
additional  deformation  which  in  turn  will  lead  to  increased  stresses  in  the  supports 
resisting  such  deformation. 

d.  By  blocking,  a  prestress  may  be  introduced  to  support  rock.  Wedging  of 
the  blocks  will  give  equal  and  opposite  forces  acting  on  the  tunnel  surface  and  the 
support  structure. 

In  the  work  reported  herein,  the  influence  of  various  material  properties  upon 
support  stresses  was  studied  assuming  load  development  primarily  through  mechanism 
(a)  described  above.  The  system  in  Figure  VI- 1  (case  a)  was  analyzed  to  study  the 
effect  of  variation  of  parameters  and  to  rank  them  in  order  of  importance.  Calculations 


were  made  varying  eaeh  oi  the  parameters  between  the  limits  given  one  at  a  time. 

The  'constant'  values  of  other  parameters  were  taken  to  be  the  midpoint  of  the  speci¬ 
fied  range.  The  combination  of  parameter  values  corresponding  to  the  worst  stress 
in  the  supports  was  then  used  to  analyze  eases  b,  c,  d  shown  in  Figures  VI-2  to 
Figure  VI-4. 

6.2.2.  Assumptions  Made  in  the  Analysis 

a.  Extent  of  the  Finite  Model 

When  an  underground  opening  is  excavated,  ehanges  in  the  stress  field  and 
associated  deformations  oeeur  in  thp  entire  rook  mass.  The  principle  of  loeal  action 
implies  that  these  ehanges  diminish  with  increasing  distance  from  the  opening.  In 
the  finite  element  model,  a  finite  region  is  generally  considered.  On  the  boundary 
of  this  finite  model,  force  or  displacement  boundary  conditions  have  to  be  prescribed. 
These  can  be  based  on  the  assumption  either  of  no  change  in  the  stress  field  or  of  no 
deformation.  Neither  of  the  two  is  true  for  finite  distances  from  the  opening  and  the 
two  assumptions  in  fact  give  bounds  to  the  correet  solution.  Nair  (1968)  and  Kulhawy 
(1972),  among  others,  have  studied  the  effeet  of  lateral  dimension  of  the  finite  model 
and  of  the  ehoiee  of  boundary  conditions  on  stresses  and  deformations  in  the  vicinity 
of  underground  openings.  For  the  present  study,  allowing  foi  the  specified  hydrostatie 
initial  stress  field,  a  preliminary  analysis  showed  that  it  would  be  adequate  to  model  a 
region  extending  approximately  seven  diameters  above  the  roof  of  the  opening,  five  dia¬ 
meters  horizontally  on  each  side  of  the  center-line  and  five  diameters  below  the  invert 
of  the  opening.  This  region  is  outlined  as  ABCD  in  Figure  VI-6.  The  boundary  conditions 
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for  the  finite  element  model  are  illustrated  in  Figure  VI-7.  The  overburden  of  the 
top  858  feet  was  replaced  by  an  equivalent  vertical  load.  The  stress  field  on  the 
vertical  faces  AB,  CD  was  assumed  to  be  unaffected  by  the  sequence  of  operations 
in  the  opening  and  the  vertical  displacement  of  the  horizontal  section  BC  was  set 
equal  to  zero. 

b.  Modelling  of  Support  Structure  and  the  Sequence  of  Operations 

Starting  with  the  initial  stress  state,  the  sequence  of  excavation,  installation 
of  supports  and  load  development  was  simulated.  The  steel  supports  were  assumed 
to  be  in  plane  stress  whereas  the  rock  and  shoterete  were  in  plane  strain.  The 
spacing  of  steel  supports  was  specified  as  three  feet.  Thus  a  three  foot  length  ot 
the  tunnel  was  supported  by  each  support  ring.  Figures  VI-8  and  VI-9  show  a  typi¬ 
cal  eross-seetion  used  in  the  finite  clement  model.  The  cross-section  of  the  steel 
rib  was  represented  by  five  finite  elements  to  obtain  reasonably  good  distribution  of 
stresses  over  the  eross-seetion.  The  shoterete  and  the  rib  were  assumed  to  be 
bonded.  If  there  is  no  bond  between  shoterete  and  the  rib,  there  is  no  load  trans¬ 
ferred  through  shear  and  the  load  transfer  is  entirely  radial.  This  situation  would 
be  similar  to  ease  (d)  except  that  the  blocking  would  be  continuous.  For  eases  b,  e, 
d,  the  wooden  blocks  were  assumed  to  be  axial  members  not  transmitting  any  bending 
moments. 

The  shoterete  ,  the  steel  rib  and  the  timber  blocks  were  assigned  the  properties 
given  in  Table  VI- 1.  It  was  assumed  that  C,  <p  specified  for  the  rock  were  obtained 
in  triaxial  tests  on  eylinderieal  specimens.  It  was  assumed  that  the  rock  properties 
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FIG.  VI-8.  Longitudinal  Cross-Section  of  Tunnel  and  Representative  Slice 


do  not  change  as  a  result  of  excavation  of  the  opening.  Additional  studies  considering 
the  effect  of  such  changes  are  reported  in  section  6.2.3. 

Tn  calculating  load  acting  on  the  supports,  it  was  assumed  that,  upon  excavation, 
*he  gravity  load  due  to  the  area  shown  shaded  in  Figure  VI-  10  was  mobilized  only  fifty 
percent.  The  rest  being  supported  hy  the  unexeavated  rock  ahead  of  the  working  face. 
The  balance  fifty  percent  was  assumed  to  become  active  upon  loss  of  support  due  to 
continued  excavation,  'lids  corresponds  to  mechanism  (a)  described  in  seetion  6.2.1. 
The  geometry  of  the  shaded  area  corresponded  to  the  overbreak  line  in  Figure  VI-2, 
ease  b.  It  was  assumed  that  the  difference  in  exeavatior  profiles  in  Figure  VI-2  and 
VI- 1  represented  the  rock  load  which  could  develop  for  eases  a  and  d.  For  case  1), 
the  rock  exerting  the  gravity  load  is  shown  shaded  in  Figure  VI- 1 1,  the  extent  was 
arbitrarily  taken  as  an  average  of  3.5  feet  thickness  beyond  the  excavation  line. 
Similarly  for  ease  e,  the  extent  (Figure  VI- 12)  was  arbitrarily  taken  as  an  average 
of  3.5  feet  thickness  beyond  the  excavation  line. 

6.2.3.  Results 
a.  Preliminary 

In  studying  the  importance  of  material  parameters,  the  quantities  of  interest 
were  the  maximum  stresses  in  the  steel  member.  It  is  customary  to  study  the  axial 
and  shearing  forces  and  bending  moments  in  such  structural  components.  Hence, 
these  quantities  were  worked  out.  It  should  be  noted  that  the  basic  output  from  the 
finite  element  analysis  is  the  components  of  stress  evaluated  at  the  center  of  each  of 
the  five  elements  into  which  the  steel  member  is  divided.  The  moments  and  forces 
were  obtained  by  numerical  integration  of  the  stress  values. 

137 


Tunnel  Cross-Section  Showing  Rock  Load  for  Case  c 


Figures  VI-13  and  VI- 14  show  the  eross-secdons  at  whieh  the  stresses  as  well 
as  the  forees  and  moments  in  the  steel  rib  were  computed. 

Hock  load  development  was  assumed  to  follow  meehanism  (ai  of  seetion  6.2.  1 . 
Additional  studies  using  meehanism  (e)  and  (d)  are  eovered  in  seetion  6.2.4. 

b.  Influence  of  Material  Properties 

Four  material  parameters,  viz.,  Young's  modulus,  Poisson's  ratio,  cohesion 
and  angle  of  internal  friction  were  to  be  considered  to  establish  their  order  ot  import¬ 
ance  in  terms  of  their  influence  on  the  support  forees. 

For  the  specified  range  uf  the  roek  properties,  the  eohes'  n  and  the  angle  of  in¬ 
ternal  friction  did  not  influence  the  stresses  in  the  structural  supports.  Figure  VI- 15 


shows  a  plot  of  Jj  -  o-^  +  ^  +  tr^,  the  first  invariant  of  the  stress  tensor  as  the  abeissa 

T  2  _  \  <ffl  -  *2' 2  +  <*2  -  0-3)2  +  (o-3  -  vtf  \  2  tl 
and  Jo  -  < - >  ,  the  seeond  invariant  of  the  stress 

I  6  \ 

deviation  tensor  as  the  ordinate.  The  generalized  Mohr-Coulomb  yield  criterion  corre¬ 
sponding  to  the  oreseribed  range  of  values  of  cohesion  and  angle  of  internal  frietion  are  shown 
and  also  the  stress  paths  traced  by  points  of  eritieal  locations  around  the  tunnel  opening. 

The  generalized  Mohr-Coulomb  yield  law  is, 

j 

a  Ji  +  J2~  k 


where: 


2  sin  i 


X/3-  (3  -  r7n0) 


k  =  6  C 

"V3“(3 


sin^) 


<P  =  angle  of  internal  frietion 
C  =  cohesion. 
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C  and  <P  are  determined  from  a  triaxial  test  on  a  cylindrical  specimen  under 
axisymmetric  radial  stress.  Derivation  of  the  above  equations  is  given  by  Singh  (1972) 

In  Figure  VI-15,  paths  A,B,C,D,  K  are  traced  by  elements  located  around  the 
face  of  the  underground  opening.  Points  A,C  refer  to  elements  at  the  invert  and  the 
crown,  respectively,  and  points  B,  D,  E  correspond  to  elements  on  the  side  of  the  open 
ing.  The  locations  are  indicated  in  Figure  VI- 13.  The  initial  state  in  all  eases  is  of 
hydrostatic  stress.  The  terminal  points  represent  the  state  after  excavation.  The 
development  of  roek  load  has  little  influence  on  the  stress  state  in  rock  for  the  speci¬ 
fied  values  of  Young's  modulus  and  Poisson's  ratio.  For  all  locations,  the  entire 
stress  history  was  found  to  be  well  below  the  yield  criterion. 

Figures  VI- 16  through  VI- 18  show  the  influence  of  variation  of  Young's  modu¬ 
lus  anon  bending  moments  and  upon  axial  and  shearing  forces.  Only  sections 
with  the  worst  forces  have  been  plotted.  The  maximum  bending  moments  were 
at  section  16  (crown),  the  maximum  axial  force  at  section  5  (side)  and  the  maximum 
shearing  force  at  section  10.  In  all  eases,  decrease  in  elastic  modulus  of  roek  was 
seen  to  result  in  increased  moments  and  forces.  This  is  indeed  to  be  expected.  The 
load  transfer  to  the  structural  support  is  dependent  upon  the  tendency  of  the  rock  face 
to  deform.  Higher  modulus  implies  less  deformation  for  the  same  roek  load  and 
consequently  less  load  transferred  to  the  steel  member.  The  decrease  in  support 
forces, with  increasing  Young's  modulus, is  rapid  at  first  and  then  is  less  pronounced. 
This  is  due  to  the  fact  that  roek  deformation  is  proportional  to  reciprocal  of  the  modu¬ 
lus.  Also,  for  very  large  moduli,  the  strains  are  extremely  small  and  difficulties 
arise  with  eemputational  precision. 
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Figures  VI-19  through  VI-211  show  the  influence  of  the  variation  in  the  Poisson's 
ratio  on  the  moments  and  axial  and  shear  forees  in  the  supports.  Again,  only 
results  for  the  worst  sections  have  been  plotted.  In  the  erown  and  the  side  (section 
10  and  5),  an  increase  in  Poisson's  ratio  results  in  decreased  bending  moments 
whereas  at  section  10  the  bending  moment  increases  somewhat.  The  higher  Poisson 
ratio  is  associated  with  a  redistribution  of  stress  in  rock.  This  redistribution  is 
reflected  in  a  more  uniform  stress  distribution  in  the  steel  support. 

For  the  specified  range  of  values  for  the  material  parameters,  the  bending 
moments  and  the  shearing  forees  were  insignificant.  The  major  effect  was  the  axial 
force  in  the  member.  Figure  VI-22  shows  the  distribution  of  longitudinal  stress  for 
the  mean  values  of  rock  parameters.  An  explanation  for  the  bending  effect  being  very 
small  may  be  found  in  the  asssumption  that  the  steel  member  is  restrained  by  the  shot- 
erete.  Significant  bending  of  the  steel  member  must  involve  significant  changes  in 
curvature.  This  is  not  possible  for  a  member  restrained  from  radial  movement  by 
relatively  unyielding  rock.  Also,  the  roek  properties  are  assumed  to  be  unaffected 
by  the  excavation  process.  This  may  not  be  true.  There  ean  be  considerable  change 
in  deformability  of  a  roek  mans  due  to  excavation  of  the  underground  opening. 

It  is  not  possible  to  assign  ranks  to  the  parameters.  It  is  clear  however,  that 
supports  in  lower  elastic  modulus  roek  will  develop  greater  stresses.  High  value  of 
Poisson's  ratio  has  the  effect  of  redistributing  stresses;  decreasing  the  peaks  and 


increasing  the  lowest  values. 


Influence  of  Poisson'G  Ratio  on  Bending  Moments 
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L  sing  the  lowest  value,  in  the  prescribed  range,  of  Young's  modulus  and  the 
minimum  as  well  as  the  maximum  values  of  Poisson's  ratio,  stresses  for  eases  a,  b, 
e,  <1  (Figures  VI-  1  through  VI-4)  of  different  excavation  profiles  and  blocking  arrange¬ 
ments  were  determined.  Figures  VI-23  through  VI-2G  show  the  distribution  of  longi- 
tundial  stresses  on  critical  scetion.  A  study  was  also  carried  out  for  ease  (a)  in 
which  no  load  transfer  through  shear  between  shoterete  and  the  rib  was  allowed  (ribs 
unbonded  to  the  shoterete).  Table  \  1-2  shows  a  comparison  of  forces  for  the  eases 
(a)  and  (d).  It  was  seen  that  the  bending  moments  in  case  (d)  were  greater  than  those 
in  ease  (a)  in  all  portions  of  the  rib  except  the  invert  section.  Also  the  axial  forces 
in  these  two  eases  are  almost  the  same  at  the  invert  section  but  were  different  at 
other  sections.  These  differences  are  due  to  the  restraint  offered  by  the  shoterete 
to  radial  deformation  of  the  rib  and  due  to  the  transfer  of  stresses  between  the  rib 
and  the  shoterete  through  shear.  Results  for  the  ease  of  no  load  transfer  through 
shear  between  the  shoterete  and  rib  showed  that  the  effect  of  bonding  of  rib  to  shot¬ 
erete  on  forces  in  the  invert  section  is  insignificant.  The  axial  forces  in  the  other 
parts  of  the  steel  rib  were  greater  for  the  case  where  load  transfer  through  bond 
was  permitted.  As  would  be  expected,  in  ease  (e),  where  the  tunnel  cross- section 
and  the  rock  load  weie  unsymmetrieal,  the  forces  in  the  steel  rib  were  also  unsym- 
metrieal.  The  axial  force  had  a  minimum  value  at  section  22  and  increased  towards 
the  invert  on  both  sides  of  section  22  to  1137  lbs.  at  seetion  1  and  913  lbs.  at  scetion 
31.  (Figure  VI-14  gives  location  of  these  sections  of  the  steel  ribs). 

6.2.4.  Additional  Studies 

Additional  studies  using  values  of  Young's  modulus  less  than  the  minimum  of  the 
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range  given  in  Table  VI-1  were  carried  out  to  consider  situations  where  the  elastic 
modulus  of  roek  in  the  vieinity  of  an  underground  opening  may  be  significantly  .educed 
by  damage  during  exeavation  or  deterioration  with  exposure  over  a  long  time.  The 
simulation  assumed  that  after  installation  » f  supports  and  placement  of  shoterete, 

Young's  modulus  may  reduce  from  an  initial  value  of  1  x  10®  pounds  per  square  inch. 

C  A 

Different  terminal  values  of  Young's  modulus  used  were  0.75  x  10  ,  0.67  x  10  ,  0.4  x 
10®,  0.25  x  10®  pounds  per  square  inch.  Poisson's  ratio  was  assumed  to  be  0.49 
throughout.  For  case  (a),  the  results  are  plotted  in  Figures  VI-27  through  VI-29. 

As  might  be  expected,  greater  reduction  in  the  clastic  modulus  was  associated  with 
greater  bending  moments  and  axial  and  shear  forces.  For  reduction  of  Young's  modu¬ 
lus  to  0.25  x  10®  psi,  the  maximum  longitudinal  stress  would  be  over  20,000  psi  (Figure 
VI-30).  Cons' dering  that  deterioration  of  roek  is  more  likely  to  occur  in  eases  b  and  d, 
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analysis  were  performed  corresponding  to  a  reduction  in  Young's  modulus  from  1x10 
to  0.4  x  10®  pounds  per  squa  *e  inch.  Figure  VI-31  shows  the  distribution  of  longitudinal 

stresses  at  critical  sections  for  the  cases  a,  b  and  d. 

Another  study  considered  mechanism  (d)  (section  6.2.1),  In  this,  using  ease  b 

of  Figure  VI-2,  the  loner  struts  or  blocking  points  were  assumed  to  be  one-dimensional 
elements  under  a  eompressive  stress  of  200  pounds  per  square  inch.  Young's  modulus  and 
Poisson's  ratio  of  rock  were  taken  to  be  1  x  10®  psi  and  0.49  respectively.  The  struct¬ 
ural  steel  member  was  determined  to  h^ve  maximum  axial  force  of  25,520  pounds  at 
section  12,  maximum  shearing  foree  of  1,604  pounds  at  section  8  and  maximum  bend¬ 
ing  moment  of  28,700  pounds  ineh  at  section  4.  The  maximum  longitudinal  stress  was 
2,750  psi  in  an  element  loeated  in  the  flange  of  section  4.  Distribution  of  longitudinal 
and  shearing  stresses  at  critical  sections  is  shown  in  Figure  VI-32. 
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-30.  Stress  Distribution  at  Critical  Sections  Case  a  (Reduction  of  E  from  1.0  x  10  to  0.25  x  10  psi) 


2750  psi 


EXTERIOR 


LONGITUDINAL  STRESS 


SHEAR  STRESS 


SECTION  4 


SECTION  8 


FIG.  VI-  32.  Distribution  of  Longitudinal  and  Shearing  Stresses  on  Critical  Sections 

Case  (I))  -  200  psi  Stress  in  Struts 


CHAPTER  VII  -  DISCUSSION 


Two  <1  life* re  nt  types  of  mechanical  behavior  of  rock  were  considered  One 
treats  rock  as  an  isotropic  elastic-plastic  material  following  a  generalized  Mohr- 
Coulomb  law  and  in  the  other,  the  rock  is  an  isotropic  linear  elastic  brittle  material 
subject  to  fracture  in  accordance  with  Griffith's  theory  or  the  Modified  Griffith 
theory.  These  two  types  of  behavior  are  representative  of  a  wide  class  of  roek 
materials. 

for  analysis  ol  stresses,  delorniation  and  progressive  failure  of  nonhomo- 
genous  fissured  rock,  the  finite  element  method  is  the  most  suitable.  In  the  pre¬ 
sent  research  effort,  previous  attempts  at  finite  element  modeling  of  the  two  types 
of  material  behavior  (elastic-plastic  and  elastied>rittle)  were  critically  examined 
and  their  limitations  noted.  The  present  development  eliminates  several  of  the 
limitations.  Analysis  of  elastic-plastic-  roek  is  based  on  Felippa's  formulation  of 
the  ineremental  stress/incrcmcntal  strain  relationship.  This  is  more  general  than 
and  includes  Prager-Naghdi's,  Reyes'  and  Reyes-Dee  re's,  and  Yamada's  formula¬ 
tions.  In  numerical  procedures,  the  incremental  approach  was  favored  because  of 
the  poor  convergence  characteristics  of  the  so  called  initial  strain  and  initial  stress 
methods.  The  procedures  developed  allow  for  incremental  conslruction/exeavation, 
arbitrary  initial  stresses,  arbitrary  geometry  and  considerable  nonhomo  ^cne  it  \  of 
of  material.  Nonmonotonic  loading  associated  with  sequential  excavation  was 


properly  allowed  for. 

For  analysis  o  progressive  failure  of  elastic-brittle  rock,  Griffith's  theory 
and  the  Modified  Grif  it*’  hcory  were  used.  Previous  work  used  the  no  tension 
approach  which  is  incorrect  and  unrealistic.  In  the  present  work,  propagation  ol 
fracture  was  considered  as  a  sequential  phenomenon.  The  stress-redistribution 
associated  with  crack  extension  was  allowed  for  as  incremental  cracking  occured 
Preexisting  joints, whether  open  or  closed, were  considered.  In  all  previous  develop¬ 
ment,  the  location  and  orientation  of  fracture  has  to  be  known  for  a  study  of  its  pro¬ 
pagation.  In  the  present  development,  fractures  can  initiate  at  randomly  oriented 
Griffith  flaws  assumed  to  be  pre-exlsteni  everywhere  in  rock.  Propagation 
follows  initiation  depending  upon  the  Griffith  criteria  being  satisfied. 

The  mathematical  models  of  material  behavior  were  fitted  into  appropriate 
variational  formulations  of  the  incremental  clastostatic  problem  for  the  case  of 
progressive  fracture  of  clastic-brittle  rock  and  the  incremental  plastieitj  problem 
for  the  case  of  progressive  failure  of  elastic-plastic  rock.  Discretization  of  the 
governing  functional  by  the  finite  element  method  gave  the  set  of  matrix  equations 
leading  to  the  problem  solution.  The  procedural  details  of  the  finite  clement  method 
are  well  known  and  therefore,  were  not  reproduced  in  the  present  report. 

The  procedures  developed  were  verified  against  existing  theoretical  solu¬ 
tions  and  experimental  data.  Excellent  agreement  was  observed.  Some  typical 
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problems  in  rock  mechanics  and  other  areas  were  solved  as  llustrative  applica¬ 
tions. 

The  techniques  developed  are  applicable  to  analysis  of  progressive  failure 
through  plastic  yield  or  brittle  fractuie.  In  either  case,  arbitrary  initial  stress 
state,  arbitrary  geometrical  configuration,  arbitrary  sequence  of  excavation/con¬ 
struction,  can  be  considered  and  a  history  of  sequential  failure  or  fracture  obtained. 
The  methods  can  allow  for  pre-existing  joints  and  fissures  and  are  applicable  to 
comparative  stability  studies  based  on  stresses  and  deformations  associated  with 
excavation  operations,  evaluation  of  structural  supports  and  loads  on  underground 
supports,  safety  analysis  of  openings,  study  of  blasting  effectiveness  under  certain 
conditions,  evaluation  of  alternative  mining  sequences  to  obtain  the  safest  construc¬ 
tion  sequence,  etc. 

The  development  was  applied  to  a  parametric  study  of  the  influence  of  rock 

properties  on  the  stresses  in  steel  supports  using  data  supplied  by  the  sp-  nsor. 

The  experimental  phase  of  the  research  program  was  concerned  with  develop¬ 
ment  of  modeling  material  so  tnat  materials  with  predetermined  shape  of  the  stress- 
strain  curve  could  be  produced  in  the  laboratory.  It  should  now  be  possible  to  carry 
model  tests  on  simulated  rock  to  predict  actual  behavior  at  site  as  also  to  verify 
computational  procedures.  Certain  assumptions  made  in  the  theory  of  plasticity 
arc  somewhat  arbitrary.  Model  experiments  would  serve  to  verify  them  and  permit 
evaluation  of  their  influence  on  the  significant  parameters  affecting  stability  of  under¬ 


ground  openings. 
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APPENDIX  A.  GRIFFITH  THEORY  OF  BRITTLE  FRACTURE 


A.  1 .  The  Original  Griffith  Theory  of  Brittle  Fracture 

I'wo  different  failure  criteria  were  proposed  by  Griffith.  They  arc  the  energy 
criterion  and  the  stress  criterion. 

A.  1. 1 .  The  Energy  Criterion 

An  extended  minimum  potential  energy  theorem  can  be  stated  as: 

If  an  clastic  solid  body  is  deformed  by  certain  specified  boundary 
forces,  the  sum  of  the  potential  energy  of  the  applied  forces  and  the  strain 
energy  stored  in  the  body  will  be  either  decreased  or  unaltered  by  the  crea¬ 
tion  of  cracks  with  traction-free  surfaces. 

Based  upon  this  theorem  and  assuming  the  creation  of  new  surfaces  required 

only  certain  amount  of  surface  energy,  the  failure  criterion  can  be  written  as: 

The  surface  energy  increased  must  be  equal  to  the  strain  energy 
released  such  that  the  foregoing  theorem  is  satisifed. 

Mathematically  the  statement  can  be  expressed  by  the  following  equation: 

-jr  <w  -  s>  -  o  (A.  n 

where  e  half  length  of  the  new  crack 

W  =  strain  energy  released 

S  -  surface  energy  required  to  form  the  new  crack. 

Eq.  (A.l)  is  considered  as  a  necessary  condition  for  stable  fracture  propagation. 

It  is  understandable  that  stress  is  highly  concentrated  at  crack  tip  when  the 
radius  of  curvature  of  the  crack  tip  is  very  small.  The  stress  concentration  may  be 
counted  for  the  failure  of  material  at  certain  stress  level  lower  than  the  strength  ci 

A-l 

- - - - - - -  .  ..  . .  . . 


the  material.  One  major  drawback  of  this  failure  criterion  is  that  the  critical  orienta¬ 
tion  of  the  Initial  flaw  which  happens  to  be  one  of  the  important  roles  in  the  fracturing 
process,  can  not  be  generally  defined.  Following  the  energy  criterion,  the  crack  would 
propagate  in  the  direction  of  the  initial  flaw.  Apparently,  this  is  just  a  special  ease  of 
fhe  general  fracturing  mode. 


A. 1.2.  Stress  Criterion 


The  stress  criterion  can  be  derived  from  the  solution  for  stress  distribution 
around  the  erack.  Following  Inglis  (1913),  the  tangential  stress  on  the  boundary  of 
an  elliptic  erack  is  given  by  (Fig.  A.l) 


/  2  ^ 

°tj  |  (W  sinh  2^o  +  (°1  -tr2)  [  e  °  cos  2  (P  V)  -  eos  20]}/ 


(eosh  2  £q-  eos  2t)  ) 


(A.  2) 


where  £  ,  V  -  orthogonal  curvilinear  coordinates 
£0  £  on  the  crack  boundary 

(Tj,  0^2  -  major  and  minor  principal  stresses,  positive  in  tension 


P  =  angle  between  the  major  axis  of  the  ellipse  and  a?  axis 


Eq.  f A. 2)  has  been  obtained  on  the  assumption  that  the  crack  is  very  flat  such  that 
£  ' :  1  exists. 

From  Eq.  (A. 2),  it  can  be  shown  (Jaeger  and  Cook,  1969)  that  the  maximum 
tensile  stress  at  fee  craek  tip  is  approximated  by 


max 


[o-i  eos2  p  -  (r ^  sin2  p  -  (o-j*  cos2  p  .  a*  sin2  p)'  ]  (A  3) 
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The  condition  ^max /dp  =  0  yields 

cos  2P  -  -  (o-j  -  cr2)  /  (o-j  »  <r2)  (A>4; 

which  defines  P  issociated  with  <rmax,  provided  0j  /  «r2  . 

Eq.  (A. 4)  also  implies 

T,  •  3,2  i  0  (A.5) 

Substituting  (A.  4)  in  (A.  3)  wc  obtain 

°max  “(<ri~(r2)  /4  (!rl  1  <r2  *  (A. 6) 

Let  <rt  denote  the  uniaxial  tensile  strength  of  the  material,  and  consider  two 
extremal  cases; 


^  °7nax  2(ri/ £Qwhen  p  is  *~et  equal  to  zero  in  Ea.(A.3)  and 
the  negative  value  of  the  square  root  is  used  to  obtain  the 
largest  value  of  ormax.  is  tensile  since  4  3<r2  ^0 

andor1  >  an  Then  (r  c.  2  a  2ir  oru  -it  n 

1  ^  max  fo  l  ^°t  t  l  u* 


(ii) 


(T 

max 


2  a  2 


when  p  — 
2 


(A. 7) 


Mathematically  (r2  may  not  be  tensile  since  Eq.  (A.5)  is  satisfied  alright,  but  physically, 
failure  could  not  happen  under  this  single  condition.  Therefore  this  becomes  a  limiting 
condition  for  Eq.  (A. 6).  It  then  follows  from  Eq.  (A. 6)  that 

2<rt  =  ((ri  -v2  / 4  <<v  v 

°r  "V1/-  8  \  <<rl  +  V  '  0  (A.  8) 

Eqs.  (A. 4),  (A.5),  (A. 7)  and  (A. 8)  constitute  the  failure  criterion.  They  are 


A -3 


summarized  as  follows; 


a-!  -  (T 2  >0 


(i)  3  <r  i  i  <r2  >  0 

<rt  -o-j  1  0 

P  "  1  (A. 9) 


OP 


3  iTj  i  cr  g  <0 


°"t  ~  J  n  -  0 


1  2 

where  crn  -  —  (0'1  -  cr2)  '  /  (cr j  »  rr2)  >0 

8 

cos  2p  -  j  (0'1  -  cr2)  /  (<r1  t  (r2) 

A .  2 .  Mocii fication  of  Griffith  Theory 

The  original  Griffith  theory  did  not  allow  for  closing  of  crack  subject  to 
compressive  stresses.  When  a  crack  closes,  sliding  may  occur  along  the  crack 
surfaces.  A  modification  permitting  crack  closure  was  proposed  by  Brace  (1960) 
and  McClintock  and  Walsh  (1962)  and  was  further  developed  by  Murrell  (1964)  and 
Hoek  and  Bieniawski  (1965). 

The  modified  theory  is  based  on  the  assumption  that  the  crack  would  close 
when  the  stress  normal  to  the  plane  of  crack  is  in  compression  and  exceeds  some 
critical  compressive  stress.  After  crack-clo  <ure,  stresses  are  carried  over  the 
crack  plane  through  interface  friction.  The  material  strength  is  increased  and 
stress  concentration  at  crack  tip  is  reduced. 


A-4 


The  conditions  are  stated  as 


"n  =  0  if  Kf  4  Kc| 


and 


'  I  °y  ~  °c  I  if  ry  >  "e 


(A.  10) 


and 

where 


The  interface  friction  and  the  effective  interface  frietion  are  given  by 

Tf  *  -'“'n  “  I  "y-'cl  •  lf  %  >ao 

T  '  Txy  -  Tf  Txy-'‘l°'y  -  Jcl 
if  ~  stress  normal  to  the  crack 


n 


(/> .  1 1  a) 
(A.  lib) 


effective  normal  stress 
cre  critical  i  >rmal  stress 

p  coefficient  of  frietion 

-  shearing  stress  along  the  eraek 

Following  the  similar  procedure  for  deriving  the  stress  criterion  and  making 


use  of  the  condition 


OT? 


=  0  ,  we  find 


2<cc  (o-  7  <  C  '  r,2) 


and 


^max  4-  ^e  4  4  T*  > 

^o 


T  =  20-  +  («re  /  <rt  +  1)  2 


2  <r* 


4 


(A.  12) 
(A .  1 3a) 

(A.  13b) 


Eq.  (A.  lib)  ean  be  rewritten  as 


2  2 

t  -  *  (Ci  -  cr2)  sin  23  +  p  (cr ^  eos  3  4  <r2  sin  3  -  <r  ) 


(A.  14) 


The  condition  t/^  3  -  0  gives 


tan  23  -  -  —  ,  provided  o-,  /  vr 


(A.  15) 
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Substituting  Eq.  (A.  15)  in  Eq.  (A.  14) 

,  2J 

T  ^  -2  “  tr2)  <  1  ^  P  M  (ffj  +  ^  "  2(rc 


Eliminating  r  between  Eq.  (A.  13b)  and  Eq.  (A.  16) 


i 

.2,2 


(OH]  -  <r2)  (  1  4  n  )  ♦  \i  (<T]  f  cr2)  -  4  o-t  (--  +  l)'2  t-  2p 


(A.  17) 


If  «rc  is  negligible,  Eq.  (A.  17)  reduces  to 

1 

(OH]  -  or2)  (  1  >  Ji2  )2+  /i  (o-j  +  <r2)  =  4  trt  (A.  18) 

Eq.  (A.]8)  coincides  with  Muhr-Coulomb  criterion  if  2  crj.  =T,  where  T  is  the  co¬ 
hesion  of  the  material.  Eq.  (A.  17)  indicates  that  the  modified  criterion  is  of  three- 
parameter  type  and  Eq.  (A.  18)  gives  a  linear  relationship  between  the  principal 
stresses. 

Another  form  of  modified  criterion  was  suggested  by  Hook  and  Bieniawski  (1966) 
and  Bieniawski  (1967): 

°2  =  trl  [  +  f  M)  /  (  VT ~2  -  n)  ]  -  tr'c  (A.  19) 

where  «r'c  -  uniaxial  compressive  strength  of  the  material.  Eq.  (A.  19)  has  the 


advantage  of  using  (r'c  which  is  much  easier  to  measure  than  tr* 


4 PPENDK  B.  development  of  test  material 

B.  1 .  The  Experimental  Program 

The  experimental  program  war  intended  to  design  laboratory  materials  whieh 
would  simulate  the  stress-strain  behavior  of  a  wide  variety  of  rock  types.  Considera¬ 
tion  would  be  given  to  brittle,  piastie  and  strain  hardening  response  both  undralned 
and  drained  pore  water  conditions.  Since  the  theoretical  solutions  would  be  expected 
to  operate  with  any  consistent  level  of  engineering  parameters  for  both  models  and 
prototypes,  it  was  not  necessary  to  develop  a  true  rock-llte  material.  Instead,  soils 
stabilized  with  additives  such  as  hydrated  lime  could  be  used.  Materials  of  this  type 
may  be  constructed  so  as  to  show  failure  response  ranging  from  clastic  brittle  to 
elastie-worlt  hardening  plastic.  The  variation  may  come  about  as  a  function  of  the 
consolidation  pressure  or  of  additive  concentration.  A  special  reason  for  usir-  soil 
to  model  rock  behavior  war  that  the  testing  and  modelling  could  be  done  using  standard 

soil  mechanics  laboratory  equipment  with  considerable  saving  in  equipment  costs  and 
development  time. 

B.2.  The  Study  Material 

The  study  material  wasabrown  silt  from  the  lacustrine  deposits  south  of  Cleveland, 
Ohio.  The  oven-dried  soil  was  combined  with  either  hydrated  lime  or  Portland  cement. 
After  adding  water,  the  soil  was  mixed  and  placed  in  a  vaeuum  extruder.  The  auger 
expelled  the  soil  in  a  saturated  condition  through  a  final  die,  ready  shaped  for  testing. 

Samples  were  sealed  in  four  layers  of  piastie  and  wax  and  then  cured  a  humid  room 
for  period  of  4  to  G  weeks. 


B-l 


A  scries  of  consolidated-draincd  triaxial  tests  indicated  that  a  4^  lime-silt 


mixture  might  have  the  desired  properties. 


As  shown  in  Figure  B-I,  the  full  range 


of  stress- strain  response  from 


elastic-brittle  to  strain  Hardening  is  available  as  a 


function  of  consolidation  pressure.  The  stress  levels  ore  within  the  usable  range  of 
the  plane  strain  device  or  of  modelling  applications.  The  tenting  program  was  lam¬ 
inated  before  plane  strain  tests  could  be  conducted.  The  different  boundary  conditions 
in  that  test  might  have  required  modifications  in  the  mix. 


B.  3.  The  Plane  Strain  Device 

A  plane  strain  device  was  purchased  from  the  Massachusetts  Institute  of 
Technology  soils  moratory.  It  has  been  shown  to  be  satisfactory  for  the  purpose  of 
the  tests  intended.  A  rectangular  sample  3.5  inch  X  3.5  inch  x  1 .4  inch  is  restrained 
so  as  to  allow  the  application  of  vert, cal  stress  by  means  of  a  piston,  control  of  stress 
or  strain  in  an  orthogonal  direction  and  the  measurement  of  stress  in  the  third  ortho¬ 
gonal  direction  against  a  fixed  plate.  The  system  is  designed  to  operate  under  a  hydro 
lie  back-pressure  and  measure  pore-pressures  during  the  test.  Schematic  drawings  in 
Figures  B-II  to  B-IV  give  details  of  the  device  and  show  elements  of  the  sample  loading 
and  arrangement.  A  complete  description  of  the  deviee  and  its  use  are  given  by  Bovee 
and  Ladd  (19701.  Hue  to  a  curtailment  of  the  scope  of  work,  no  tests  were  carried  out 


under  the  current  research  program. 
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FIG.  B-I.  Stress-Strain  Behavior  of  4%  Line  f’.iit  (T ri axial  Consolidated 
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FIXED  END  PLATENS 
(SIDE  VIEW  OF  DEVICE) 

FIG.  B-II 


(From  Bovee-1970) 
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HORIZONTAL  LOADING  SYSTEM 
(END  VIEW  OF  DEVICE) 


(From  Bovec- 1970) 
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